footnote

# A graph counterexample to davies’ conjecture

Gady Kozma

###### Abstract.

There exists a graph with two vertices $x$ and $y$ such that the ratio of the
heat kernels $p(x,x;t)/p(y,y;t)$ does not converge as $t\to\infty$.

This paper is concerned with a conjecture of Brian Davies from 1997 on the
heat kernel of Riemannian manifolds [D97, §5]. We will not disprove the
conjecture as stated, but rather transform it to the realm of graphs using a
well-known (though informal) “dictionary” between these two categories, and
build a graph that will serve as a counterexample. We will make some remarks
on how the construction might be carried over back to the category of
manifolds in the end, but we will not give all details. The bulk of this paper
is about graphs.

We start by describing the conjecture in its original setting. Let $M$ be a
connected Riemannian manifold, and let $p$ be the _heat kernel_ associated
with the _Laplace-Beltrami operator_ on $M$. Then Davies conjecture states
that for any $M$ and any 3 points $x,y,z\in M$ the limit

$\lim_{t\to\infty}\frac{p(x,y;t)}{p(z,z;t)}$ (1)

exists and is positive. Here $p(x,y;t)$ is the value of the heat kernel at
time $t$ and at points $x$ and $y$. This property is known as the “strong
ratio limit property” (where the “weak” version is an averaged result due to
Döblin, [D38]) or SRLP for short. So Davies’ conjecture is that in these
settings SRLP always holds. SRLP holds for manifolds with one end [D82] and
for strongly Liouville manifolds (i.e. manifolds where any positive harmonic
function is constant), see [ABJ02, corollary 2.7] who also make interesting
connections between these properties and the infinite Brownian loop.

Ratio limit properties were considered for Markov chains even earlier. If $M$
is any Markov chain on a countable state space, then we say that $M$ satisfies
SRLP if (1) holds for any three states $x$, $y$ and $z$, where $p(x,y;t)$ is
the probability that the Markov chain started at $x$ will be at $y$ at time
$t$. For general Markov chains there are a few examples where SRLP does not
hold. Clearly it does not hold when the Markov chain has some kind of
periodicity. F. J. Dyson constructed an example of an aperiodic recurrent
Markov chain which does not satisfy SRLP [C60, part I, §10]. That example
utilizes long chains of states with only one outgoing edge, which the walker
must traverse sequentially. In particular it is not _reversible_.

Now, the Laplace-Beltrami operator is self-adjoint so a proper analog of
Davies conjecture needs to assume that the Markov chain is reversible.
Reversible Markov chains are also known as random walks on weighted graphs.
The issue of periodicity can be dealt with by looking at random walk in
continuous time or at lazy random walk. Lazy random walk is a walk where the
walker, at every step, chooses with probability $\frac{1}{2}$ to stay where it
is, and with probability $\frac{1}{2}$ moves to one of the neighbours (with
probability proportional to the weights).

The result is this paper is

######  Theorem.

There exists a connected graph $G$ with bounded weights and two vertices
$x,y\in G$ such that the heat kernel of the lazy random walk satisfies

$\frac{p(x,x;t)}{p(y,y;t)}\nrightarrow$ (2)

as $t\to\infty$.

Let us remark on the “bounded weights” clause. When doing analogies between
manifolds and graphs, it is often assumed that the manifold has bounded
geometry and the graph has bounded weights. Davies, however, explicitly does
not make the assumption of bounded geometry. Thus one might wonder what is
exactly the graph analog. All this is mute, of course, since the
counterexample does has bounded weights (and hence a manifold example
constructed along the same line should have bounded geometry).

It is easy to see that in these setting (reversible, irreducible, lazy) this
ratio must be bounded between two constants independent of $n$. Hence if it
does not converge then it must fluctuate between two values. The proof
constructs the graph with two halves, denoted $H^{e}$ and $H^{o}$ ($e$ and $o$
standing for even and odd, $H$ standing for half), which are connected by one
edge, $(x,y)$. On the “odd scales”, $H^{e}$ will “look like $\mathbb{Z}^{22}$”
while $H^{o}$ will “look like $\mathbb{Z}^{3}$”. This means that to get from
$x$ to $x$ (where $x$ is on the $H^{e}$ side), the most beneficial strategy is
to move to $y$ as fast as possible, spend most of your time on the $H^{o}$
side and return to the $H^{e}$ side only in the last minute. Clearly this
would mean that $p(x,x;t)$ is smaller than $p(y,y;t)$ as the random walk
starting from $y$ can stay on its side at all time, not losing the constant
that $x$ needs for the maneuver. On the “even scales” the picture is reversed
and $y$ is at a disadvantage. See figure 1 — drawing in 22 dimensions might
have distracted the reader, so the figure demonstrates the construction in 1
and 2 dimensions. The smallest two braces in the figure are the first scale,
in which $H^{o}$ is really one dimensional and $H^{e}$ is really two-
dimensional. The larger two braces indicate the second scale. This time
$H^{o}$ is a network of lines so it should be thought about as two
dimensional, while $H^{e}$ is a thick column, so it should be thought about as
one dimensional. The third scale is only hinted at in this figure, but one can
imagine that $H^{o}$ now becomes a thick band, so it is again one-dimensional,
while $H^{e}$ becomes a network of these thick columns and bands, so it is
back to being two dimensional.

$x$$y$$H^{e}$$H^{o}$

Figure 1. The graphs $H^{o}$ and $H^{e}$

As one might expect, the numbers $3$ and $22$ have no particular significance.
They both have to be $>2$, since otherwise our graph would be recurrent and
recurrent graphs always satisfy SRLP [O61, theorem 3]. And of course they have
to be different. We took here the large value $22$ in order to be able to be
wasteful in various points (sum over times and such stuff), but the proof
could proceed with any value larger than $3$.

This paper was first written in 2006. I wish to take this opportunity to
apologize to all those who has to wait so long for it to appear, with no real
reason. My intentions were good but my time management was abysmal. I wish to
thank Yehuda Pinchover for telling me about the problem and for reading early
drafts. Partially supported by the Israel Science Foundation.

##  1\. Proof

The construction uses $\mathbb{Z}^{d}$-like graphs as building blocks, so we
start with quoting a few results on these. We first recall the notion of rough
isometry [K85].

######  Definition.

Let $X$ and $Y$ be two metric spaces. We say that $X$ and $Y$ are roughly
isometric if there is a constant $C$ and a map $\varphi:X\to Y$ with the
following properties:

  1. (i)

For all $x$ and $y$ in $X$,

$\frac{1}{C}d(x,y)-C\leq d(\varphi(x),\varphi(y))\leq Cd(x,y)+C$

  2. (ii)

The image of $\varphi$ is roughly dense, i.e. for all $y\in Y$ there is an
$x\in X$ such that $d(y,\varphi(x))\leq C$

If $G$ and $H$ are graphs we say that they are roughly isometric if they are
roughly isometric when considered with the metric $d$ being the graph
distance, namely $d(x,y)$ is the length of the shortest path between $x$ and
$y$, or $\infty$ if no such path exists.

With this definition we can state the following standard result, essentially
due to Delmotte.

######  Lemma 1.

Let $G$ be a graph roughly isometric to $\mathbb{Z}^{d}$. Then the heat kernel
$p$ for the lazy walk on $G$ satisfies, for all $t\geq 1$,

$ct^{-d/2}\leq p(x,x;t)\leq Ct^{-d/2}.$

Here $G$ is a simple graph — we do not allow weights or multiple edges. $C$
and $c$ are constant which do not depend on $t$. In general we will use $c$
for constants which are small enough and $C$ for constants which are large
enough, and different appearance of $c$ and $C$ might relate to different
constants.

###### Proof.

By Delmotte’s theorem [D99] any $G$ which satisfies volume-doubling and the
Poincaré inequality, satisfies

$p(x,x;t)\approx\frac{1}{|B(x,\sqrt{t})|}$

where $B(x,r)$ is the ball around $x$ with radius $r$ (again with the graph
distance), and $|B(x,r)|$ is the sum of the degrees of the vertices in $B$.
The notation $X\approx Y$ is short for $cY\leq X\leq CY$. The fact that $G$ is
roughly-isometric to $\mathbb{Z}^{d}$ gives

$|B(x,r)|\approx r^{d}$ (3)

so we would get $p(x,x;t)\approx t^{-d/2}$, as needed. So we need only show
that $G$ satisfies volume doubling and Poincaré inequality.

Now, the definition of volume doubling is that for every $x$ a vertex of $G$
and every $r\geq 1$

$|B(x,2r)|\leq C|B(x,r)|$

and it follows immediately from (3). The Poincaré inequality is not much more
complicated. The definition is: for every vertex $x$, for every $r$ and for
every function $f:B(x,2r)\to\mathbb{R}$,

$\sum_{y\in B(x,r)}\deg(y)|f(y)-\overline{f}|^{2}\leq Cr^{2}\sum_{(y,z)\in
E(B(x,2r))}(f(y)-f(z))^{2}$ (4)

where

$\overline{f}=\frac{1}{|B(x,r)|}\sum_{y\in B(x,r)}\deg(y)f(y).$

and where $\deg(y)$ is the degree of $y$, and $E(B)$ is the set of edges both
whose vertices are in $B$. Now, $\mathbb{Z}^{d}$ satisfies the Poincaré
inequality (see e.g. [PSC, §4.1.1]). It is well-known and not difficult to see
that the Poincaré inequality is preserved by rough isometries (it uses the
fact that $\sum\deg(y)|f(y)-a|^{2}$ is minimized when $a=\overline{f}$). This
finishes the proof.∎

######  Lemma 2.

Let $G$ be a graph roughly-isometric to $\mathbb{Z}^{d}$, $d\geq 3$ and let
$x$ be some vertex. Let $p$ be the probability that lazy random walk starting
from $x$ returns to $x$ for the first time at time $t$. Then $p\geq
ct^{-d/2}$.

###### Proof.

Let $p_{1}$ be the same probability but without the restriction that this is
the first return to $x$. This is exactly $p(x,x;t)$ and by theorem 1 we have
$p_{1}\geq ct^{-d/2}$. Fix some $K$ and examine the event that the random walk
returns to $x$ at $t$ and also at some time $s\in[K,t-K]$. Let $p_{2}$ be its
probability. Using the other direction in theorem 1 we can write

$p_{2}\leq\sum_{s=K}^{t-K}p(x,x;s)p(x,x;t-s)\leq
C\sum_{s=K}^{t-K}s^{-d/2}(t-s)^{-d/2}\leq CK^{1-d/2}t^{-d/2}.$

Since $d\geq 3$ we can choose $K$ sufficiently large such that
$p_{2}\leq\frac{1}{2}p_{1}$. So we know that with probability $p_{1}-p_{2}\geq
ct^{-d/2}$ the walk does not return to $x$ between $K$ and $t-K$. If it does
reach $x$ before time $K$, do some local modification so that it does not. For
example, if the original walker reached $x$ at some time $s<K$ and on the next
step went to some neighbour $y$ of $x$, modify it to walk to $y$ in the first
step and stay there for $s$ steps (remember that our walk is lazy) and then
continue like the original walker. Clearly this “costs” only a constant and
ensures our walker does not visit $x$ in the interval $[1,K]$. Do the same for
the interval $[t-K,t-1]$, losing another constant. The finishes the lemma.∎

######  Lemma 3.

Let $G$ be a graph and let $p$ be the heat kernel for the lazy walk on $G$.
Let $t$ and $s$ satisfy that $|t-s|\leq\sqrt{t}$. Then

$|p(x,x;t)-p(x,x;s)|\leq
C\frac{|t-s|\log^{3}t}{\sqrt{t}}p(x,x;t)+Ce^{-c\log^{2}t}.$

###### Proof.

Denote by $q(x,y;t)$ the heat kernel for the _simple_ random walk on $G$. Then
by definition

$p(x,x;t)=\sum_{i=0}^{t}q(x,x;i)\binom{t}{i}2^{-t}.$

Writing the same formula for $p(x,x;s)$ and subtracting we get

$|p(x,x;t)-p(x,x;s)|\leq\sum_{i}q(x,x;i)\left(\binom{t}{i}2^{-t}-\binom{s}{i}2^{-s}\right)=\Sigma_{1}+\Sigma_{2}$

where $\Sigma_{1}$ is the sum over all $|t-2i|\leq\sqrt{t}\log t$ and
$\Sigma_{2}$ is the rest. A simple calculation with Stirling’s formula shows
that

$2^{-t}\binom{t}{i}=\sqrt{\frac{2}{\pi
t}}\exp\left(-\frac{(t-2i)^{2}}{2t}\left(1+O\left(\frac{|t-2i|+1}{t}\right)\right)\right)$

And with some more calculations

$\Sigma_{1}\leq C\sum_{|t-2i|\leq\sqrt{t}\log t}q(x,x;i)\sqrt{\frac{2}{\pi
t}}e^{-(t-2i)^{2}/2t}\frac{|t-s|\log^{3}t}{\sqrt{t}}\leq
C\frac{|t-s|\log^{3}t}{\sqrt{t}}p(x,x;t)$

while

$\Sigma_{2}\leq C\sum_{|t-2i|>\sqrt{t}\log t}e^{-c(t-2i)^{2}/t}\leq
Ce^{-c\log^{2}t}$

proving the lemma. ∎

###### Proof of the theorem.

Abusing notations, for subsets $H\subset\mathbb{Z}^{d}$ we will not
distinguish between $H$ as a set and as an induced subgraph of
$\mathbb{Z}^{d}$ ($d$ will be $22$). For the construction we need a
sufficiently fast increasing sequence $a_{1}<a_{2}<\dotsb$. We further assume
that $a_{k}$ are all even and that $a_{k-1}$ divides $\frac{1}{2}a_{k}$. It
would have probably been enough to choose $a_{k}=2^{a_{k-1}}$, but it turns
out simpler to choose the $a_{k}$ inductively, and we perform this as follows.
Let $a_{1}=2$. Assume now $a_{1},\dotsc,a_{k-1}$ have been defined. Define,
for integers $m<\frac{1}{2}l$ and $i\in\\{1,\dotsc,22\\}$,

$\displaystyle Q_{l,m,i}$
$\displaystyle:=\left\\{\vec{n}\in\mathbb{Z}^{22}:\left|n_{i}\textrm{ mod
}l\right|\leq m\right\\}$ $\displaystyle Q_{l,m}$
$\displaystyle:=\bigcup_{\begin{subarray}{c}I\subset\\{1,\dotsc,22\\}\\\
|I|=19\end{subarray}}\bigcap_{i\in I}Q_{l,m,i}.$

Here $n\textrm{ mod
}l\in\\{-\lfloor\frac{l-1}{2}\rfloor,\dotsc,\lfloor\frac{l}{2}\rfloor\\}$. In
words, $Q_{l,m,i}$ is a $21$-dimensional subspace of $\mathbb{Z}^{22}$
orthogonal to one of the axes, fattened up by $2m+1$ (a “slab”) and repeated
periodically with period $l$. $Q_{l,m}$ is the collection of all
$3$-dimensional subspaces, fattened and repeated similarly. The particular
point $\vec{0}$ is in fact contained in all $\binom{22}{3}$ of these
$3$-dimensional slabs which will be a little inconvenient, so let us shift
$Q_{l,m}$ by

$v(m)=\Big{(}\underbrace{{\textstyle\frac{1}{2}}m,\dotsc,{\textstyle\frac{1}{2}}m}_{3\textrm{
times}},\underbrace{\vphantom{{\textstyle\frac{1}{2}}m}0,\dotsc,0}_{19\textrm{
times}}\Big{)}.$

In the shifted set $Q_{l,m}+v(m)$ the geometry of the neighbourhood of
$\vec{0}$ is simpler, it is contained in just one slab. Compare to the figure
on page 1. The point $x$ is in the _middle_ of a fat column and not at the
intersection of a column and a band.

We want to use these graphs with $l=a_{j}$ and $m$ a little larger than
$a_{j-1}$. Precisely, define

$b_{j}=\sum_{k=1}^{j}a_{k}.$

With this choice of $b_{j}$, $Q(a_{j},b_{j-1},i)$ contains only complete
components of $Q(a_{l},\linebreak[4]b_{l-1},i)$ for each $l<j$. Each such
component is either contained in $Q(a_{j},b_{j-1},i)$ or disjoint from it. The
same holds for the translations $Q(a_{l},b_{l-1},i)+v(a_{l})$ (we need here
that $a_{l}>4a_{l-1}$ so let us assume this from now on). For brevity, define
$v_{j}=v(a_{j})$.

We may now define two graphs, denoted by $H_{k-1}^{\textrm{e}}$ and
$H_{k-1}^{o}$ (“e” and “o” standing for even and odd) by

$H_{k-1}^{\textrm{e/o}}:=\bigcap_{\begin{subarray}{c}2\leq j\leq k-1\\\
j\textrm{ even/odd}\end{subarray}}(Q_{a_{j},b_{j-1}}+v_{j}).$

We shall usually suppress the $k-1$ from the notation. It is not difficult to
check that $H^{\textrm{e/o}}$ are both roughly isometric to $\mathbb{Z}^{22}$
(the rough isometry constant depends on the “past” $a_{1},\dotsc,a_{k-1})$.
Therefore by lemma 1 we see that there exists an $\alpha$ (again, depending on
the past) such that

$p_{H^{\textrm{e/o}}}(x,x;t)\leq\alpha t^{-11}.$ (5)

Examine now the graphs

$F_{k-1}^{\textrm{e/o}}:=H_{k-1}^{\textrm{e/o}}\cap\left\\{\vec{n}\in\mathbb{Z}^{22}:\left|n_{i}\right|\leq
b_{k-1}\>\forall i=4,\dotsc,22\right\\}.$

$F^{\textrm{e/o}}$ are both roughly isometric to $\mathbb{Z}^{3}$ so by lemma
2 there exists some $\beta$ such that

$\mathbb{P}_{F^{\textrm{e/o}}}(\mbox{the walk returns to }\vec{0}\mbox{ for
the first time at }t)\geq\frac{1}{\beta}t^{-3/2}.$ (6)

Define $\gamma_{k}:=\left\lceil\max\\{\alpha,\beta\\}\right\rceil$ (as usual,
$\left\lceil\cdot\right\rceil$ stands for the upper integer value). With these
we can define $a_{k}$ to be any even number satisfying
$a_{k}>2\gamma_{k}^{4}+4a_{k-1}$ and such that $a_{k-1}$ divides
$\frac{1}{2}a_{k}$. This completes the description of the induction, and we
define

$H_{\infty}^{\textrm{e/o}}:=\bigcap_{\begin{subarray}{c}2\leq j\\\ j\textrm{
even/odd}\end{subarray}}(Q_{a_{j},a_{j-1}}+v_{j}).$

These graphs will be the two halves of our target graph $G$.

Before continuing, let us collect some simple facts about
$H_{\infty}^{\textrm{e/o}}$:

  1. (i)

$H_{\infty}^{\textrm{e/o}}$ is connected — in fact we used this indirectly
when we claimed $H_{k}^{\textrm{e/o}}$ are roughly isometric to
$\mathbb{Z}^{22}$.

  2. (ii)

$H_{\infty}^{\textrm{e/o}}$ are transient — this follows because each contains
a copy of $\mathbb{Z}^{3}$ (namely $\\{n_{4}=\dotsb=n_{22}=0\\}$) and
transience is preserved on adding edges. This last fact follows from
conductance arguments, see e.g. [DS84].

Define therefore the escape probabilities

$\varepsilon^{\textrm{e/o}}:=\mathbb{P}_{H_{\infty}^{\textrm{e/o}}}^{\vec{0}}(R(t)\neq\vec{0}\,\forall
t>0)$

($R$ being the random walk on the graph) and let
$\delta:=\frac{1}{2}\min\\{\varepsilon^{\textrm{e}},\varepsilon^{\textrm{o}}\\}$.
Define the graph $G$ by connecting $H_{\infty}^{\textrm{e}}$ to
$H_{\infty}^{\textrm{o}}$ with a single edge between the two $\vec{0}$ with
weight $\delta$. Define $x:=\vec{0}^{\textrm{e}}$ and
$y=\vec{0}^{\textrm{o}}$. This is our construction and we need to show (2),
which will follow if we show that, for $k$ sufficiently large,

$\left.\begin{aligned} p(x,x;t_{2k})&\geq 3p(y,y;t_{2k})\\\
p(x,x;t_{2k+1})&\leq{\textstyle\frac{1}{3}}p(y,y;t_{2k+1})\end{aligned}\right\\}\quad
t_{k}:=\gamma_{k}^{4}.$ (7)

We will only prove the even case, the odd will follow similarly.

Examine therefore $p(x,x;t_{2k})$. Since $a_{2k}>t_{2k}$ we get that

$H_{\infty}^{\textrm{e/o}}\cap[-t_{2k},t_{2k}]^{22}=H_{2k}^{\textrm{e/o}}\cap[-t_{2k},t_{2k}]^{22}$

or in other words, the steps after $2k$ do not effect us at all. Similarly it
is possible to simplify the last stage namely

$\displaystyle H_{2k}^{\textrm{e}}\cap[-t_{2k},t_{2k}]^{22}$
$\displaystyle=H_{2k-1}^{\textrm{e}}\cap(Q_{a_{2k},b_{2k-1}}+v_{2k})\cap[-t_{2k},t_{2k}]^{22}=$
$\displaystyle=H_{2k-1}^{\textrm{e}}\cap\left\\{\vec{n}\in\mathbb{Z}^{22}:\left|n_{i}\right|\leq
b_{2k-1}\>\forall i=4,\dotsc,22\right\\}=F_{2k-1}^{\textrm{e}}$

(here is where these translations by $v_{j}$ are used). By (6),

$\displaystyle p_{G}(x,x;t_{2k})$
$\displaystyle\geq\frac{1}{2}\mathbb{P}_{H_{2k}^{e}}(R\mbox{ returns to
}x\mbox{ for the first time at }t)\geq$
$\displaystyle\stackrel{{\scriptstyle\textrm{(\ref{eq:beta})}}}{{\geq}}\frac{1}{2\gamma_{2k}}t_{2k}^{-3/2}=\frac{1}{2}t_{2k}^{-7/4}$
(8)

(the $\frac{1}{2}$ comes from the first step).

To estimate $p(y,y;t_{2k})$ we divide the event $\\{R(t_{2k})=y\\}$ according
to whether $R$ “essentially goes through $x$” or not. Formally, denote by
$T_{1}$ and $T_{2}$ the first and last time before $t_{2k}$ when $R(T)=x$ (if
this does not happen, denote $T_{1}=\infty$ and $T_{2}=-\infty$). Then we
define

$\displaystyle p_{1}:=\mathbb{P}(T_{1}>\gamma_{2k},\,R(t_{2k})=y)\qquad
p_{2}:=\mathbb{P}(T_{2}<t_{2k}-\gamma_{2k},\,R(t_{2k})=y)$ $\displaystyle
p_{3}:=\mathbb{P}(T_{1}\leq\gamma_{2k},\,T_{2}\geq
t_{2k}-\gamma_{2k},\,R(t_{2k})=y)$

so that $p(y,y;t_{2k})\leq p_{1}+p_{2}+p_{3}$.

Now, $p_{1}$ and $p_{2}$ are easy to estimate. As above we have

$H_{2k}^{\textrm{o}}\cap[-t_{2k},t_{2k}]^{22}=H_{2k-1}^{\textrm{o}}\cap[-t_{2k},t_{2k}]^{22}$

so (5) applies and we get

$\mathbb{P}_{H_{\infty}^{\textrm{o}}}^{y}(R(t)=y)\leq\gamma_{2k}t^{-11}\quad\forall
t\leq t_{2k}.$ (9)

Therefore

$\displaystyle p_{1}$
$\displaystyle\leq\sum_{t=\gamma_{2k}}^{t_{2k}-1}\mathbb{P}(T_{1}=t,\,R(t_{2k})=y)+\mathbb{P}(T_{1}=\infty,\,R(t_{2k})=y)\leq$
$\displaystyle\leq\sum_{t=\gamma_{2k}-1}^{t_{2k}-2}\mathbb{P}_{H_{\infty}^{\textrm{o}}}(R(t)=y)+\mathbb{P}_{H_{\infty}^{\textrm{o}}}(R(t_{2k})=y)\leq$
$\displaystyle\stackrel{{\scriptstyle\textrm{(\ref{eq:yy})}}}{{\leq}}\sum_{t=\gamma_{2k}-1}^{t_{2k}-2}\gamma_{2k}\cdot
t^{-11}+\gamma_{2k}\cdot t_{2k}^{-11}\leq
C\gamma_{2k}^{-9}=Ct_{2k}^{-9/4}\stackrel{{\scriptstyle\textrm{(\ref{eq:pxxle})}}}{{=}}o(p(x,x;t))$
(10)

and similarly for $p_{2}$. As for $p_{3}$, we have

$\mathbb{P}(T_{1}\leq\gamma_{2k})\leq\Big{(}\sum_{i=0}^{\infty}\mathbb{P}_{H_{\infty}^{\textrm{o}}}(r\textrm{
visits }y\textrm{ }i\textrm{ times before
}\gamma_{2k})\Big{)}\cdot\delta\leq\frac{\delta}{\epsilon^{\textrm{o}}}\leq\frac{1}{2}$

and similarly (using time reversal) for $\mathbb{P}(T_{2}\geq
t_{2k}-\gamma_{2k})$. Hence we get

$p_{3}\leq\frac{1}{4}\max_{t_{2k}-2\gamma_{2k}\leq s\leq t_{2k}}p(x,x;s)$

and by lemma 3,

$\displaystyle p_{3}$
$\displaystyle\leq\frac{1}{4}p(x,x;t_{2k})\left(1+O\left(\frac{\gamma_{2k}\log^{3}t_{2k}}{\sqrt{t_{2k}}}\right)\right)+O(e^{-c\log^{2}t_{2k}})\leq$
$\displaystyle\stackrel{{\scriptstyle\textrm{(\ref{eq:pxxle})}}}{{\leq}}\frac{1}{4}p(x,x;t_{2k})(1+o(1)).$

With the estimate (10) for $p_{1}$ and the corresponding estimate for $p_{2}$
we get

$p(y,y;t_{2k})\leq p(x,x;t_{2k})\left(\frac{1}{4}+o(1)\right).$

A completely symmetric argument shows that at $t_{2k+1}$ the opposite occurs:

$p(x,x;t_{2k+1})\leq p(y,y;t_{2k+1})\left(\frac{1}{4}+o(1)\right)$

proving the theorem.∎

######  Remark.

If you want an example with unweighted graphs, this is not a problem. $H^{e}$
and $H^{o}$ are already unweighted, so the only thing needed is to connect
them, instead of with an edge of weight $\delta$, with a segment sufficiently
long such that the probability to traverse it is $\leq\delta$. The proof
remains essentially the same.

##  2\. Manifolds

We would like to show a Manifold $M$ and two points $x,y\in M$ such that the
heat kernel on $M$ satisfies

$\frac{p(x,x;t)}{p(y,y;t)}\nrightarrow$

as $t\to\infty$. Here is how one might translate the construction of our
theorem to the settings of manifolds. The dimension of the manifold plays
little role, so we might as well construct a surface.

For a subset $H\subset\mathbb{Z}^{22}$ one can associate a manifold $H^{*}$ by
replacing each vertex $v\in H$ with a sphere $v^{*}$ and every edge with a
empty, baseless cylinder. Since the degree of every vertex in $H$ is $\leq 44$
we may simply designate 44 disjoint circles on $\mathbb{S}^{2}$ and attach the
cylinders to the spheres at these circles. This is reminiscent of the well-
known “infinite jungle gym” construction, see some lovely pictures in [ON].
The exact method of doing so is unimportant since anyway the manifold that we
get is roughly isometric to $H$, considered as an induced subgraph of
$\mathbb{Z}^{22}$ (one of the nice features of rough isometry is that
continuous and discrete objects may be roughly isometric, rough isometry
inspects only the large scale geometry). Clearly $H^{*}$ can be made
$C^{\infty}$.

One can then construct a (possibly different) sequence $a_{k}$ and two
manifolds $\big{(}H_{\infty}^{\textrm{e/o}}\big{)}^{*}$ with the only
difference is that the $\alpha$ and $\beta$ must satisfy (5) and (6) for our
choice of the ∗ operation. This should be possible since
$\big{(}H_{k}^{\textrm{e/o}}\big{)}^{*}$ and
$\big{(}F_{k}^{\textrm{e/o}}\big{)}^{*}$ are roughly isometric to
$\mathbb{Z}^{22}$ and $\mathbb{Z}^{3}$ respectively. Instead of Delmotte one
can use the manifold version [SC95] (or rather, [D99] is the graph version of
earlier results for manifolds, see [SC95] for historical remarks).

The argument for the transience of
$\big{(}H_{\infty}^{\textrm{e/o}}\big{)}^{*}$ should also be direct
translation. Each contains a submanifold (with boundary) which is roughly
isometric to $\mathbb{Z}^{3}$ and therefore is transient. Since transience is
equivalent to the fact that for some $c>0$ every function which is $1$ at $x$
and $0$ at infinity satisfies that the Dirichlet form $\langle\nabla f,\nabla
f\rangle>c$, and since restricting to a submanifold only decreases the
Dirichlet form, we see that $\big{(}H_{\infty}^{\textrm{e/o}}\big{)}^{*}$ are
transient. Denote by

$\epsilon^{\textrm{e/o}}=\inf_{x\in
v^{*},v\sim\vec{0}^{\textrm{e/o}}}\mathbb{P}^{x}\left(W[0,\infty)\cap\big{(}\vec{0}^{\textrm{e/o}}\big{)}^{*}=\emptyset\right)$

where $W$ here is the Brownian motion on the manifold
$\big{(}H_{\infty}^{\textrm{e/o}}\big{)}^{*}$; and where the infimum is taken
over all $x$ belonging to a sphere $v^{*}$ where $v$ is some neighbor of
$\vec{0}^{\textrm{e/o}}$ in $H_{\infty}^{\textrm{e/o}}$. One can now define
$\delta=\frac{1}{2}\min(\varepsilon^{\textrm{e}},\varepsilon^{\textrm{o}})$
and connect $\vec{0}^{\textrm{e}}$ to $\vec{0}^{\textrm{o}}$ by a cylinder
sufficiently thin (or sufficiently long) such that the probability to traverse
it in either direction before reaching a neighboring sphere is $\leq\delta$.
This concludes a possible construction of a manifold $M$, and one may take $x$
to be an arbitrary point in $\big{(}\vec{0}^{\textrm{e}}\big{)}^{*}$ and $y$
and arbitrary point in $\big{(}\vec{0}^{\textrm{o}}\big{)}^{*}$.

The proof that $p(x,x;t)/p(y,y;t)$ does not converge should not require
significant changes. We note that in our case it is possible for a Brownian
motion at time $t$ to exit the box $[-t,t]^{22}$, but it is exponentially
difficult to do so. Hence, for example, instead of (8) we get

$p(x,x;t_{2K})\geq\frac{1}{2\gamma_{2k}^{2}}t_{2k}^{-3/2}-Ce^{-ct_{2k}}\geq\frac{1}{4}t_{2k}^{-7/4}$

for $k$ sufficiently large. Another point to note is that lemma 3 needs to be
replaced with an appropriate analog.

## References

  * [ABJ02] Jean-Philippe Anker, Philippe Bougerol and Thierry Jeulin, _The infinite Brownian loop on a symmetric space_. Rev. Mat. Iberoamericana 18:1 (2002), 41–97. Available at: projecteuclid.org
  * [C60] Kai Lai Chung, _Markov chains with stationary transition probabilities_. Die Grundlehren der mathematischen Wissenschaften 104, Springer-Verlag, 1960.
  * [D82] Brian E. Davies, _Metastable states of symmetric Markov semigroups. II_. J. London Math. Soc. 26:3 (1982), 541–556. Available at: oxfordjournals.org
  * [D97] Brian E. Davies, _Non-Gaussian aspects of heat kernel behaviour_. J. London Math. Soc. 55:1 (1997), 105–125. Available at: oxfordjournals.org
  * [D38] Wolfgang Doeblin, _Sur deux problèmes de M. Kolmogoroff concernant les chaînes dénombrables_. Bull. de la Soc. Math. de France, 66 (1938), 210–220. Available at: numdam.org
  * [DS84] Peter G. Doyle and Laurie J. Snell, _Random walks and electric networks_. Carus Mathematical Monographs, 22. Mathematical Association of America, Washington DC, 1984. Available at: dartmouth.edu/~doyle
  * [D99] Thierry Delmotte, _Parabolic Harnack inequality and estimates of Markov chains on graphs_ , Rev. Mat. Iberoamericana 15:1 (1999), 181–232. Available at: rsme.es
  * [K85] Masahiko Kanai, _Rough isometries, and combinatorial approximations of geometries of non-compact Riemannian manifolds_. J. Math. Soc. Japan 37:3 (1985), 391–413. projecteuclid.org
  * [ON] Barrett O’Neill, _An infinite jungle gym_. Available at: ucla.edu/~bon
  * [O61] Steven Orey, _Strong ratio limit property_. Bull. Amer. Math. Soc. 67:5 (1961), 571–574. Available at: ams.org
  * [PSC] Christophe Pittet and Laurent Saloff-Coste, _A survey on the relationships between volume growth, isoperimetry, and the behavior of simple random walk on Cayley graphs, with examples_. Available from: cornell.edu/~lsc
  * [SC95] Laurent Saloff-Coste, _Parabolic Harnack inequality for divergence-form second-order differential operators_. Potential Anal. 4:4 (special issue, 1995), 429–467. Available at: springerlink.com


# Pathway Histogram Analysis of Trajectories: A general strategy for
quantification of molecular mechanisms

Ernesto Suárez<EMAIL_ADDRESS>National Cancer Institute,
Frederick, Maryland  Daniel M. Zuckerman<EMAIL_ADDRESS>OHSU Center for
Spatial Systems Biomedicine and Department of Biomedical Engineering, Oregon
Health and Science University, Portland, Oregon

###### Abstract

A key overall goal of biomolecular simulations is the characterization of
“mechanism” – the pathways through configuration space of processes such as
conformational transitions and binding. Some amount of heterogeneity is
intrinsic to the ensemble of pathways, in direct analogy to thermal
configurational ensembles. Quantification of that heterogeneity is essential
to a complete understanding of mechanism. We propose a general approach for
characterizing path ensembles based on mapping individual trajectories into
pathway classes whose populations and uncertainties can be analyzed as an
ordinary histogram, providing a quantitative “fingerprint” of mechanism. In
contrast to prior flux-based analyses used for discrete-state models,
stochastic deviations from average behavior are explicitly included via direct
classification of trajectories. The histogram approach, furthermore, is
applicable to analysis of continuous trajectories. It enables straightforward
comparison between ensembles produced by different methods or under different
conditions. To implement the formulation, we develop approaches for
classifying trajectories, including a clustering-based approach suitable for
both continuous-space (e.g., molecular dynamics) or discrete-state (e.g.,
Markov state model) trajectories, as well as a “fundamental sequence” approach
tailored for discrete-state trajectories but also applicable to continuous
trajectories through a mapping process. We apply the pathway histogram
analysis to a toy model and an extremely long atomistic molecular dynamics
trajectory of protein folding.

###### keywords:

folding, kinetics, first-passage time, Markov state model, non-Markov analysis

MSM, Markov state model; NM, non-Markovian; RMSD, root mean-squared deviation

##  1 Introduction

Protein conformational processes, which include folding and conformational
transitions, are essential to innumerable cellular functions 1, and a key goal
of computer simulations is to quantify such processes 2, 3, 4. Characterizing
the “mechanisms” (pathways) of such processes, the series of conformational
states through which a system transits during the process of interest, is at
the heart of the biomolecular simulation enterprise.

Because straight-ahead molecular dynamics (MD) simulations generally extend to
timescales that are only a fraction of what is required to exhibit
biologically significant conformational processes 5, 2, specialized
computational methods have been developed to infer long-time behavior from
affordable simulations. A number of “path sampling” techniques have been
developed in recent decades which generally attempt to construct continuous or
nearly continuous trajectory ensembles based on principles of non-equilibrium
statistical mechanics 3, 6, 7, 8. Further, path optimization approaches
estimate a single dominant pathway 9, 10, and based on a set of short
trajectories, Markov state models (MSMs) use statistical inference to “stitch
together” mechanistic pathways in a reduced configurational description
consisting of discrete states 11, 12, 13, 14, 15, 16, 17.

The analysis and interpretation of pathways constitutes another challenge,
since their intrinsic structure is often hidden behind the stochastic nature
of the system and the limited sampling. Our overall strategy is to classify
pathways in a well-defined way that enables (i) statistical comparison to
another path ensemble, perhaps of unbiased reference data, (ii) ready analysis
of the dominant transition mechanisms and the degree of heterogeneity in a
path ensemble, and (iii) applicability to continuous or discrete trajectories.
We always assume that the macrostates of interest, A and B (e.g., folded and
unfolded) have already been defined but the discussion applies to arbitrary
macrostates. To classify the pathway ensembles from any approach (e.g., MSM or
MD trajectories), we require a scheme that maps trajectories to pathway
classes in a physically reasonable way.

For analysis of MSMs, a flux-based pathway classification strategy 14, 18 has
been widely used. The approach uses an intuitive mapping of flux pathways that
allows the fractional flow through each path to be computed. Fluxes, however,
represent _average_ behavior and simply following flux lines excludes
stochastic but expected switching between (non-crossing) flux lines, as
described below. Of equal importance, flux lines cannot be used to _classify_
actual trajectories because of the stochastic behavior just noted. In general,
even a sequence of MSM states generated probabilistically from the MSM
transition matrix could not be uniquely assigned to a single flux path. Of
equal importance here, MD trajectories, even if mapped to a discrete MSM
space, also could not be classified by flux-based paths.

A very simple model can illustrate the complexity of the path ensemble and
suggest strategies to simplify its analysis. Suppose we have a very coarse
partition of the configurational space with only four states; a initial state
$A$, a final state $B$, and two intermediate states $I1$ and $I2$ that connect
$A$ and $B$ as shown in Fig. 1a. Following our “chemical intuition” we can
identify four different mechanisms to go from $A$ to the final state $B$ (Fig.
1b). However, the path ensemble is much more complex, having an unlimited
number of paths (examples in Fig. 1c) . We could also identify an infinite
family of paths, for instance, sequences following the pattern
$A\rightarrow(I1\rightarrow I2\rightarrow\mathrm{A}\rightarrow
I1)_{n}\rightarrow B$, for any positive integer $n$.

Figure 1: Simple model kinetic model. (a) Kinetic model with four states, the
arrows are drawn between states kinetically connected. (b) Possible mechanisms
from $A$ to $B$. (c) Examples of additional paths with “unproductive” steps.

In order to propose a mechanism we need to reduce the dimensionality of the
problem. A simplified description could be given by how the reactive flux from
$A$ to $B$ is distributed over the network of “edges” that connect the states
14, 18. However, a reactive flux description averages out potentially
important elements of the path ensemble. In our previous example, for
instance, each edge in Fig 1a would be unidirectional in a flux picture, and
the I1-I2 edge would have zero flux in a symmetric system, thus excluding
possible paths that move transversely or backwards and failing to separate
distinct pathways that may share common segments. We want a description of the
path ensemble without discarding possible paths, including sub-optimal
“orthogonal-space” transitions that occur in realistic trajectories.

To better describe the mechanism and evaluate the performance of possible
models, we propose to map the entire path ensemble to a limited number of
classes: see Fig. 2. Each class would be a family of paths with a likelihood
easy to quantify. The probability distribution (histogram) over the classes
would act like the fingerprint of the path ensemble. Once the mapping function
is defined, it would be very easy to compare the actual MD simulation with
those generated from models and evaluate their performance; we just have to
map them all with the same function. To summarize, the classification will
allow us to compare complex ensembles via the comparison of simple histograms.

Figure 2: Pathway histogram analysis of trajectories. Trajectories are mapped
to distinct pathway classes and then a histogram is constructed. The approach
provides a quantitative “fingerprint” of mechanism and enables statistical
comparison of trajectory ensembles.

In the remainder of the paper, we demonstrate how the histogramming strategy
can be applied to both discrete-state and continuous trajectories. For
discrete-state trajectories, such as from MSMs, histograms can be constructed
via direct mapping from trajectory to pathway class or indirectly by
clustering and then classifying. We introduce a direct mapping which yields
what we term the “fundamental sequence” — roughly, the discrete-state
trajectory with loops removed. We further demonstrate that pathway classes can
be generated from arbitrary clusterings of paths via a Voronoi decomposition –
where the initial clustering can be performed based on any algorithm and
distance metric. Continuous trajectories can be classified indirectly using
the clustering-to-classes approaches or by mapping the trajectories to a
discrete basis and then applying the fundamental-sequence – or any other –
classification.

##  2 Methods

Grouping together “similar” pathways is probably the simplest way to describe
the mechanism of a biomolecular process. The goal is to expose the “structure”
of the path ensemble through a regular classification procedure. In this
section we will first show an automated classification scheme based on what we
called the _fundamental sequence_ (FS) which is designed for discrete
trajectories, e.g., the trajectories we would obtain from MSMs. Then, we will
discuss how an initial clustering in path space can be used as starting point
for classification of both discrete and continuous pathways.

###  2.1 Classification based on direct mapping to pathway classes

If every transition trajectory can be directly mapped to one of a discrete set
of pathway classes, a path histogram is created. Such a histogram can be used
for mechanistic analysis in its original form, or the class bins can be
combined using a similarity metric for a more coarse grained description.

The classes themselves can be generated in different ways. They could be based
on physical intuition – for example, by considering the sequences of certain
discrete events of interest which occur during the transitions.

Here we present an automated classification procedure for trajectories in
discrete spaces based on the fundamental sequence (FS) idea. Such discrete
trajectories might result from mapping detailed configurations to a discrete
set of configuration-space regions, or from some dynamics defined in the
discrete space, such as a Markov state model. Roughly speaking, the FS of a
discrete path is what remains after eliminating the “loops” in the path. A
loop is any sub-segment of the path where the initial and final states are the
same (do not confuse with initial and final _macro_ states). In Fig. 1c, for
instance, the segment of the path with arrows in light red are loops. Even
self transitions –when the systems stays in the same state after $\tau$– can
be considered loops of size zero. If two paths share the same FS, then they
will belong to the same class. There is not, however, a unique way to remove
the loops and we need an unambiguous definition for FS.

Imaging drawing a graph $G(\mathbf{s})$ from the _single_ discrete
trajectory/path $\mathbf{s}=(s_{1},s_{2},...,s_{m})$ from $A$ to $B$
($s_{1}\in A$, $s_{m}\in B$), where the nodes of the graph would be all the
$N$ states. The graph edges are built following one single rule: if there is
one or more _direct_ transitions from $m$ to $l$ in $\mathbf{s}$ (i.e.,
$\exists\,i:s_{i}=m\;\mathrm{and}\;s_{i+1}=l$ ), a _directed_ edge
$m\rightarrow l$ is added to $G$. At the end, some edges will have a single
direction ($m\rightarrow l$), others will be bidirectional ($m\rightleftarrows
l$) and there will also be in general isolated nodes not visited by the
particular path $\mathbf{s}$. We assume all edges have known transition
probabilities $T_{ml}$, which for molecular systems can be approximated using
MD trajectories as is done in MSM construction 11, 12, 13, 14, 15, 16, 17.

The _fundamental sequence_ (FS) of a path $\mathbf{s}$ is the connected linear
subgraph $\mathbf{q}$ of nodes in $G(\mathbf{s})$ that maximizes the product
of the transition probabilities through all possible paths $\Gamma(G)$
consistent with $G(\mathbf{s})$. That is

$\mathrm{FS}(\mathbf{s})=\operatorname*{arg\,max}_{\mathbf{q}\in\Gamma(G)}\prod_{i=1}^{|\mathbf{q}|-1}T_{q_{i},q_{i+1}}.$
(1)

where $|\mathbf{q}|$ is the number of elements in the path
$\mathbf{q}\in\Gamma(G(\mathbf{s}))$, and $|\mathbf{q}|-1$ is the number of
transitions.

An equivalent and more convenient formulation would be

$\mathrm{FS}=\operatorname*{arg\,min}_{\mathbf{q}\in\Gamma(G)}\left\\{\sum_{i=1}^{|\mathbf{q}|-1}-\log
T_{q_{i},q_{i+1}}\right\\},$ (2)

where have transformed the problem of maximizing the product of the transition
probabilities to the classical problem of finding the shortest distance in a
graph, in this case, using the pseudo distance $\delta(m,l)=-\log T_{ml}$.
There are well known strategies to solve this kind of problem; here we use
Dijkstra’s algorithm 19.

A simple example can be seen in Fig. 1, where the FSs of the paths in Fig. 1c
are the corresponding sequences in Fig. 1b. In general, any possible path from
$A$ to $B$ in that simple model will correspond to one of the four
_fundamental sequences_ drawn in Fig. 1b.

For the description of the mechanism, and in particular for equilibrium
ensembles which exhibit mechanistic symmetry for forward and reverse
directions 20, it is more natural to use a partition of path space which is
symmetric by construction. That is, for a system in equilibrium, the
distribution of the folding paths over the classes should be equal to the
distribution of the unfolding paths on the same space. Instead of maximizing
the likelihood of the sequence given the connectivity of the path, we will
maximize the likelihood of the “round trip”. The symmetrized fundamental
sequence FS∗ is therefore defined as

$\mathrm{FS}^{*}=\operatorname*{arg\,min}_{\mathbf{q}\in\Gamma(G)}\left\\{\sum_{i=1}^{|\mathbf{q}|-1}-\log(T_{q_{i},q_{i+1}}T_{q_{i+1},q_{i}})\right\\}.$
(3)

From now on, when me mention the fundamental sequence, we really mean this
symmetrized version.

Lastly we note that the prescriptions above for computing the FS are not
themselves unique and could be modified to normalize for path length and/or
employ a more explicit maximum flux definition. Regardless of the definition,
such FS classes provide a means for binning arbitrary stochastic discrete-
state trajectories and creating a histogram.

###  2.2 Classification based on clustering in path space - for discrete or
continuous trajectories

To generate a path histogram, we require true _classification_ of path space
into non-overlapping, space-tiling regions, which typically is not achieved
using clustering methods 21. However, it is straightforward to use arbitrary
clusters and an arbitrary distance metric in path space to generate Voronoi-
cell classification of the space. Most simply, a representative member of each
cluster can be taken as the reference structure for a Voronoi cell defined
using the distance metric. Although the distance metric used for clustering
and Voronoi classification need not be the same, here we will assume they are
identical.

Importantly, because clustering can be performed for either continuous or
discrete trajectories, the corresponding Voronoi classification is equally
general. We will demonstrate clustering-based classification using discrete
trajectories.

It is well known that clustering methods do not yield unique results, nor
necessarily reproducible results for a stochastic method 21, and this
indeterminacy will carry over path classification. Nevertheless, any
reasonable clustering can be expected to yield reasonable classes for
mechanistic analysis and, moreover, _any_ classification will provide a
“fingerprint” enabling quantitative comparison between different ensembles.

####  2.2.1 Clustering and classification of discrete paths based on simple
state sets

For future study, here we suggest a very simple distance metric based on state
commonality for trajectories in discrete space. In fact, this approach
discards sequence information. Before computing any distance between two
trajectories, say, $\mathbf{s}=(s_{1},s_{2},...,s_{m})$ and
$\mathbf{q}=(q_{1},q_{2},...,q_{n})$ where $s_{1},q_{1}\in A$ and
$s_{m},q_{n}\in B$, we are going to transform them to the sets
$\mathcal{S}=\bigcup^{m}_{i=1}\\{s_{i}\\}$ and
$\mathcal{Q}=\bigcup^{n}_{i=1}\\{q_{i}\\}$ respectively.

Now we define the distance $d$ between $\mathcal{S}$ and $\mathcal{Q}$ as

$d(\mathcal{S},\mathcal{Q})=|(\mathcal{S}\setminus\mathcal{Q})\cup(\mathcal{Q}\setminus\mathcal{S})|,$
(4)

where $|.|$ is the cardinality operator, $\setminus$ is the set difference
(i.e., $A\setminus B=\\{x:x\in A\;\mathrm{and}\;x\notin B\\}$). In simple
words, $d(\mathcal{S},\mathcal{Q})$ is the number of non-common elements
between the sets $\mathcal{S}$ and $\mathcal{Q}$.

The distance $d$ is also called symmetric difference between the sets
$\mathcal{S}$ and $\mathcal{Q}$ and it is a true mathematical distance (i.e.,
$d(\mathcal{S},\mathcal{Q})\geq 0$,
$d(\mathcal{S},\mathcal{Q})=d(\mathcal{Q},\mathcal{S})$,
$d(\mathcal{S},\mathcal{Q})=0\Longleftrightarrow\mathcal{S}=\mathcal{Q}$ and
$d(\mathcal{S},\mathcal{Q})\geq
d(\mathcal{S},\mathcal{R})+d(\mathcal{R},\mathcal{Q})$)[ref]. However, what
makes $d$ special is that it is a simple and very efficient metric to define
Voronoi cells.

In order to generate candidates for the Voronoi centers we first do
hierarchical clustering based on $d$ and then the centroid $\zeta_{k}$ of the
cluster $C_{k}$ is selected as shown in Eq. 5. The set $\\{\zeta_{k}\\}$,
together with $d$ will completely determine to which class any path belongs
to.

$\zeta_{k}=\operatorname*{arg\,min}_{\mathbf{s}\in C_{k}}\,\max_{\mathbf{q}\in
C_{k}}d(\mathcal{S},\mathcal{Q}).$ (5)

In this case, paths are reduced to sets and, in contrast to the FS, the
distance used is independent of the ordering of states in the path. The FS
might be more intuitive and it is probably a better tool to interpret results.
Also, discrete distances are not ideal for Voronoi-based classification since
it is more likely to to find points equidistant to multiple Voronoi centers.
Nevertheless, the clustering approach is equally valid for our goal of
defining a partition of the path space that allows “fingerprinting” and
comparison of different path ensembles.

####  2.2.2 Clustering and classification of continuous trajectories

As a final example for possible future work, we could use the Hausdorff or
Fréchet metric recently explored by Beckstein and coworkers 22 to classify an
ensemble of continuous trajectories once the Voronoi centers in path space are
defined. Although here we are recommending to use clustering as an
intermediate step to define those centers, any other reasonable approach can
be valid.

###  2.3 Defining macrostates: Kinetic clustering of microstates based on
commute times

So far we have been talking about classification or clustering of pathways.
However, implicit in pathways are the macrostates which transition pathways
connect. Here we cluster microstates in order to identify or define or
macrostates, e.g., the folded and unfolded macrostates in a protein. There is
no a unique way to do it, and that problem by itself could be the subject of
another publication/research. Nevertheless, that is not the main goal of the
present work and we are just going to propose one of the multiple ways in
which the definition of macrostates can be done.

Here we employ a systematic, but fairly simple idea method building on prior
work 11, 23. The idea is to perform hierarchical kinetic clustering of MSM
microstates, denoted as $i$ or $j$, while ensuring a separation of timescales.
We can use as a distance metric the MFPT, or more precisely, the “round-trip”
time, $t_{ij}$ = MFPT($i\to j$) + MFPT($j\to i$). A slight limitation of the
implementation suggested below is the use of (inexact) Markovian MFPTs,
because the unbiased direct estimation would be very noisy.

We employ a simple procedure: a hierarchical (or progressive) clustering based
on a cutoff $t_{cut}$ and a factor $m$ defining the targeted separation of
timescales. If the round-trip time ($t_{ij}$) is less that $t_{cut}$ then we
merge the states. Here is the proposed procedure:

  1. 1.

Compute MFPT matrix $M$ and add it to $M^{T}$ to obtain the round-trip times
$\\{t_{ij}\\}$

  2. 2.

While $\min(\\{t_{ij}\\})<t_{cut}$:

     * •

Merge the corresponding states

     * •

Recompute $\\{t_{ij}\\}$ (step 1) for merged states

  3. 3.

Increase $t_{cut}$ until the maximum remaining element of $t_{ij}$ becomes
less than a pre-defined multiple $m$ (e.g., 10) of $t_{cut}$. Use as
macrostates the two clusters from the prior step with the maximum
$t_{ij}>m\,t_{cut}$

##  3 Results

###  3.1 Classification based on the FS

In this section we demonstrate the utility of the FS-based classification for
the analysis and interpretation of complex stochastic pathways. We start with
a two-dimensional toy model where it is easy to illustrate the how this method
simplifies, in a meaningful way, the path ensemble. Then, we will follow the
same steps on a much more complex system, protein B folding and unfolding from
atomistic simulation.

####  3.1.1 Toy Model

Our toy model is a two-dimensional system with an energy function given by:

$E_{2D}/k_{B}T=\left\\{\begin{array}[]{lr}\frac{3}{2}f(x,y)+\frac{9}{10^{4}}g(x,y)^{2}g(y,x)^{2}&\mbox{if}\>\;0<x,y<6\pi\\\
\infty&\mbox{otherwise,}\end{array}\right.$ (6)

where $f(x,y)=1-\sin(x)\sin(y)$ and $g(x,y)=x-y-9\sin(y/3)$. Fig. 3 show the
shape of the potential energy function $(E_{2D}/k_{B}T)$.

Figure 3: Simple 2D toy model. The energy values are in $k_{B}T$ units.

.

This is in some way a more complex version of our first example (see Fig. 1).
Intuitively, we would have two dominant pathways or “mechanisms”. However, in
order to do the classification based on the FS, the model has to be discrete.

To analyze paths in this model, we built a (discrete) Markov model based on a
(continuous-space) Monte Carlo simulation. Simulations of around $10^{6}$
steps were done on the energy surface given in Eq. 6 using trial moves
$(\delta x_{i},\delta y_{i})$ taken randomly from a uniform distribution in
the region $[-\pi,\pi]\times[-\pi,\pi]$. The Markov states were defined by the
squares $[k\pi,k\pi+\pi)\times[m\pi,m\pi+\pi)$, where $(k,m)$ is any
combination of integers in $[0,5]$. For the path analysis we are also
considering an initial state $A=\\{(x,y)\in[0,\pi)^{2}\\}$ and a final state
$B=\\{(x,y)\in[5\pi,6\pi)^{2}\\}$ (see Fig. 3).

From now on, our analysis will focus on the discrete states of Markov model
instead of the original continuous energy surface. For that reason the quality
of the Markov model as a representation of the original model is not the main
concern. Our goal for now is the classification of discrete paths and not how
they where obtained. For visualization purpose, however, we are going to map
discrete sequences back to the original energy surface.

Starting from the state $A$ (now a discrete index), $200$ discrete
trajectories were generated until they reach $B$. The discrete trajectories
are mapped, just for visualization, to the original energy surface, in
particular to the center of each state
($[k\pi,k\pi+\pi)\times[m\pi,m\pi+\pi)$). The result, after adding noise in
order avoid perfect overlap between distinct discrete trajectory fragments, is
shown in Fig. 4. Due to the stochastic nature of the paths it is almost
impossible to see some kind of “structure” in the path ensemble, even when the
model was designed to have two dominant paths.

One the other hand, the classification based on the FS can clearly simplify
the interpretation of the stochastic behavior of our model. Fig. 5 shows the
FSs plotted on the original energy surface; being the thickness of the
individual FSs proportional to their observed probability. Notice that the FSs
inherit the symmetry of the original model but we cannot expect a perfect
symmetry since they are obtained by statistically.

Figure 4: 200 simulated paths from $A$ to $B$. The paths were generated from a
discrete Markov model that _learned_ from the original trajectory and mapped
back to the original energy surface. Noise were added just for visualization
in order to avoid perfect overlap between discrete trajectories. The energy
values are in $k_{B}T$ units.

.

Figure 5: Path representation using fundamental sequences (FSs). The 200 paths
from $A$ to $B$ shown in Fig. 4 are mapped to their fundamental sequences,
plotted over the original energy surface. The thickness of the FSs is
proportional to their probabilities. The energy is in $k_{B}T$ units. Figure
6: Toy-model pathway histogram. The plot shows the summed histogram - i.e.,
the cumulative distribution of the classes defined by fundamental sequences,
the light green region represents the 95% confidence interval. Data are
derived from the 200 paths shown in Fig. 4.

The cumulative distribution of the FSs from $A$ to $B$ is shown in Fig. 6. The
pathway classes are ordered by their preponderance, i.e., the lower the index
the more probable. In this simple case, the two dominant FSs represent about
one third of the path ensemble, and only the first $20$ pathway classes have
more than one observed event. If there were another model of the system that
we wanted to evaluate in terms of the path ensemble, we would just have to
compare its distribution of the _same_ FSs, with those shown Fig. 6.

####  3.1.2 Protein B atomistic (un)folding trajectories

Next we study an atomistic, explicitly solvated model of Protein B, a much
more complex system with $47$ residues that was simulated using regular
molecular dynamics for $104\mu$s, enough time to observe multiple
folding/unfolding events 24. As noted above, trajectories in a continuous
space can be treated by a cluster-then-classify strategy applied in path
space. Here, however, we recapitulate the strategy used in the toy model of
projecting configurations into a discrete space and using the fundamental-
sequence analysis.

We started by reducing the configurational space of the system down to 40
discrete states using pyEMMA 25, yielding the space in which FSs will be
computed. We used the same states to build a non-Markovian (NM) model 26,
which can also be described as a “history-augmented MSM.” The folded and
unfolded states where defined by a cutoff $\,t_{cut}=45$ns following the
kinetic clustering approach described above.

To illustrate the value of the pathway histogram as a yardstick, we compare
the results obtained from (a) the actual MD simulation analyzed using the
discrete MSM states mentioned above, with two NM models, (b) one using all the
history available and (c) another NM model built using only 20ns of history,
as described previously26 See Fig. 7. As we can see, the NM models reproduce
reasonably well the populations of at least the first five FS, representing
almost $90\%$ of the path ensemble. The folding/unfolding mechanisms – which
are symmetric 20 – are dominated in this classification by a single FS (index
0) with almost 70% of the statistical weight.

Figure 7: Protein B pathway histogram. The cumulative distribution (summed
histogram) of the fundamental-sequence classes is shown. The reference 95%
confidence interval of the class distribution based on MD data is shown as the
green regaion. Distributions obtained from the non-Markovian model with full
history information (solid blue) and 20ns of history (dashed blue) are also
shown. Figure 8: Primary Protein B folding pathways. Representative
intermediate structures are shown from the top two fundamental-sequence
classes of Fig. 7, along with the percent of transitions taking each pathway.

Fig. 8 is a simplified representation of the two dominant pathways where each
structure – folded, unfolded and intermediate – is a representative element of
the corresponding pathway class. The intermediate structures suggest the FS
analysis has automatically extracted truly distinct pathways. We note that
percentages in Fig. 8 include data for folding only, and hence differ slightly
(due to noise) from the data of Fig. 7.

We emphasize that other classification schemes could be used, so our findings
are not unique. However, we expect the results to be structurally reasonable
given the underlying basis in the clustered states derived from pyEMMMA 25.
Future work will investigate the degree to which path histograms can be
considered “robust” – i.e., insensitive to reasonable perturbations to the
underlying featurization. In any case, the FS analysis can bring quantitative
insight to the understanding of complex dynamical processes like protein
folding, as well as providing a standard yardstick for comparisons.

##  4 Conclusions

Although transition pathways have been analyzed in the past by a number of
quantitative approaches 14, 18, 22, no universally applicable analysis has
previously been introduced, to our knowledge. Here we have provided employed a
simple statistical structure, the histogram, along with a very general
strategy for classifying trajectories into bins – a strategy applicable to
true stochastic trajectories whether in discrete or continuous space. This
histogram “fingerprint,” although not unique, enables a standardized
comparison of mechanisms predicted by different models and/or simulation
methods.

The pathway histogram can be constructed in arbitrary ways, and we have
suggested several approaches for both discrete and continuous-space
trajectories. Notably, the automated “fundamental sequence” (FS)
classification of discrete-state trajectories removes less informative pathway
loops while retaining the potentially stochastic path “backbone” that is
unlikely to follow flux lines. Continuous-space trajectories can employ the FS
classification scheme if configurations are first mapped to discrete states,
as we have done. Further, any clustering technique can be converted to a
classification scheme using a Voronoi procedure; the only requirement is a
distance metric in path space on which to base the clustering.

The pathway histogram approach is very generic and could equally be applied to
other fields studying stochastic transitions in complex systems, including
systems biology 27, 28, 29.

We greatly appreciate constructive comments on this work from John Chodera and
Jeremy Copperman. We thank D.E. Shaw Research for providing access to the
protein folding trajectory data set. This work was partially supported by NIH
Grants R01GM115805 and P41GM103712, as well as NSF Grant MCB-1119091.

## References

  * Alberts et al. 2002 Alberts, B.; Johnson, A.; Lewis, J.; Raff, M.; Roberts, K.; Walter, P. _Molecular Biology of the Cell (4th ed.)_ ; Garland Science: New York, 2002
  * Zuckerman 2011 Zuckerman, D. M. _Annu Rev Biophys_ 2011, _40_ , 41–62
  * Chong et al. 2017 Chong, L. T.; Saglam, A. S.; Zuckerman, D. M. _Curr. Opin. Struct. Biol._ 2017, _43_ , 88–94
  * Bolhius et al. 2002 Bolhius, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. _Annu. Rev. Phys. Chem._ 2002, _53_ , 291–318
  * Karplus and McCammon 2002 Karplus, M.; McCammon, J. A. _Nat Struct Biol_ 2002, _9_ , 646–652
  * Bello-Rivas and Elber 2015 Bello-Rivas, J. M.; Elber, R. _J. Chem. Phys._ 2015, _142_ , 094102
  * Warmflash et al. 2007 Warmflash, A.; Bhimalapuram, P.; Dinner, A. R. _The Journal of chemical physics_ 2007, _127_ , 114109
  * Allen et al. 2006 Allen, R. J.; Frenkel, D.; ten Wolde, P. R. _The Journal of chemical physics_ 2006, _124_ , 024102
  * Maragliano et al. 2006 Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. _J. Chem. Phys._ 2006, _125_ , 024106
  * Henkelman et al. 2000 Henkelman, G.; Uberuaga, B. P.; Jónsson, H. _J. Chem. Phys._ 2000, _113_ , 9901
  * Chodera et al. 2007 Chodera, J. D.; Singhal, N.; Pande, V. S.; Dill, K. A.; Swope, W. C. _J. Chem. Phys._ 2007, _126_ , 155101–17
  * Noé et al. 2007 Noé, F.; Horenko, I.; Schütte, C.; Smith, J. C. _J. Chem. Phys._ 2007, _126_ , 155102
  * Bowman et al. 2009 Bowman, G. R.; Beauchamp, K. A.; Boxer, G.; Pande, V. S. _J. Chem. Phys._ 2009, _131_ , 124101
  * Noé et al. 2009 Noé, F.; Schütte, C.; Vanden-Eijnden, E.; Reich, L.; Weikl, T. R. _Proc. Natl. Acad. Sci. U. S. A._ 2009, _106_ , 19011–19016
  * Bowman et al. 2010 Bowman, G. R.; Ensign, D. L.; Pande, V. S. _J. Chem. Theory Comput._ 2010, _6_ , 787–794
  * Prinz et al. 2011 Prinz, J.-H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; Schütte, C.; Noé, F. _J. Chem. Phys._ 2011, _134_ , 174105
  * Bowman et al. 2014 Bowman, G. R.; Pande, V. S.; Noé, F. _An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation_ ; Springer: Dordrecht, Netherlands, 2014
  * Berezhkovskii et al. 2009 Berezhkovskii, A.; Hummer, G.; Szabo, A. _The Journal of Chemical Physics_ 2009, _130_
  * Dijkstra 1959 Dijkstra, E. W. _Numerische Mathematik_ 1959, _1_ , 269–271
  * Bhatt and Zuckerman 2011 Bhatt, D.; Zuckerman, D. M. _Journal of chemical theory and computation_ 2011, _7_ , 2520–2527
  * Barbakh et al. 2009 Barbakh, W. A.; Wu, Y.; Fyfe, C. _Non-Standard Parameter Adaptation for Exploratory Data Analysis_ ; Springer, 2009; pp 7–28
  * Seyler et al. 2015 Seyler, S. L.; Kumar, A.; Thorpe, M. F.; Beckstein, O. _PLOS Computational Biology_ 2015, _11_ , 1–37
  * Zhang et al. 2010 Zhang, X.; Bhatt, D.; Zuckerman, D. M. _Journal of chemical theory and computation_ 2010, _6_ , 3048–3057
  * Lindorff-Larsen et al. 2011 Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. _Science_ 2011, _334_ , 517–520
  * Scherer et al. 2015 Scherer, M. K.; Trendelkamp-Schroer, B.; Paul, F.; Pérez-Hernández, G.; Hoffmann, M.; Plattner, N.; Wehmeyer, C.; Prinz, J.-H.; Noé, F. _Journal of chemical theory and computation_ 2015, _11_ , 5525–5542
  * Suárez et al. 2016 Suárez, E.; Adelman, J. L.; Zuckerman, D. M. _Journal of chemical theory and computation_ 2016, _12_ , 3473–3481
  * Morelli et al. 2009 Morelli, M. J.; ten Wolde, P. R.; Allen, R. J. _Proceedings of the National Academy of Sciences_ 2009, _106_ , 8101–8106
  * Donovan et al. 2016 Donovan, R. M.; Tapia, J.-J.; Sullivan, D. P.; Faeder, J. R.; Murphy, R. F.; Dittrich, M.; Zuckerman, D. M. _PLoS computational biology_ 2016, _12_ , e1004611
  * Tse et al. 2018 Tse, J., Margaret; Chu, B. K.; Gallivan, C. P.; Read, E. L. _PLoS computational biology_ 2018, _14_ , e1006336


# Studying the bound state of the $K^{-}p$ system in the Bethe-Salpeter
formalism

Chao Wang<EMAIL_ADDRESS>Center for Ecological and Environmental
Sciences, Key Laboratory for Space Bioscience $\&$ Biotechnology, Northwestern
Polytechnical University, Xi’an 710072, China  Liang-Liang Liu
<EMAIL_ADDRESS>College of Physics and Information Engineering,
Shanxi Normal University, Linfen 041004, China  College of Nuclear Science and
Technology, Beijing Normal University, Beijing 100875, China  Xin-Heng Guo
<EMAIL_ADDRESS>College of Nuclear Science and Technology, Beijing Normal
University, Beijing 100875, China

###### Abstract

We study the $s$-wave kaon-nucleon bound state with the strangeness $S=-1$ in
the Bethe-Salpeter formalism in the ladder and instantaneous approximations.
We solve the Bethe-Salpeter equation of the bound state and obtain the Bethe-
Salpeter amplitude. It is shown that the $K^{-}p$ bound state exists in this
formalism. We also study the decay width of the bound state based on the
Bethe-Salpeter techniques. The mass of this bound state is 1422 MeV and its
decay width is obviously smaller than that of $\Lambda(1405)$. These results
indicate that there may be some other structures in the observed resonance.

###### pacs:

11.30.Er, 12.39.-x, 13.25.Hw

##  I Introduction

In the last decade, several nonconventional states were observed K.A. ;
Guo:2017jvc . In the light baryon spectrum, one of such resonances is
$\Lambda$(1405), which is just below the $\bar{K}N$ threshold and emerges in
the meson-baryon scattering amplitude with $I(J^{P})=0(1/2^{-})$ and
strangeness $S=-1$ K.A. ; Hyodo:2011ur . It was discovered in the $\Sigma\pi$
invariant mass spectrum of the channel $K^{-}p\to\pi\pi\pi\Sigma$
Dalitz:1959dn ; Dalitz:1960du ; Alston:1961zzd . Various experimental and
theoretical investigations for $\Lambda(1405)$ have been preformed in recent
years. The failure of the traditional quark model to interpret $\Lambda(1405)$
makes people consider exotic configurations such as the $\bar{K}N$ molecular
structure Dalitz:1967fp ; Bohm:1977mz . This molecular configuration has been
supported by recent lattice QCD results Hall:2014uca .

Meanwhile, several studies propose that $\Lambda$(1405) has a two-pole
structure and the spectrum in experiments exhibits one effective resonance
shape Jido:2003cb . Recently, the SIDDHARTA Collaboration has determined the
energy shift and width of kaonic hydrogen, which provides a strong and direct
constraint on the $\bar{K}N$ scattering amplitude Bazzi:2011zj ; Bazzi:2012eq
. Based on this constraint, various studies confirm the two-pole structure of
$\Lambda$(1405) Ikeda:2011pi ; Ikeda:2012au ; Mai:2012dt ; Guo:2012vv ;
Mai:2014xna . In this work, we will focus on the bound state below the
$\bar{K}N$ threshold which might consist of the proton and antikaon.

The chiral perturbation theory (ChPT) is a useful effective field theory to
deal with meson-meson (baryon) interactions at the low energy scale
Gasser:1984gg ; Ecker:1994gg . A systematic and successful approach (the
unitary chiral approach) combining ChPT and the unitarity condition of the
scattering amplitude has been developed to describe the $K^{-}p$ scattering
data Hyodo:2011ur ; Oset:1997it . In this paper, we will consider the Bethe-
Salpeter equation method. The unitary chiral approach is based on the
scattering theory. In this approach, the scattering amplitude can be obtained
by regarding the two-body interaction obtained in chiral perturbation theory
as the potential and solving the Lippmann-Schwinger equation Hyodo:2011ur ;
Oset:1997it . The two-particle Bethe-Salpeter equation is derived from the
relativistic quantum field theory. The basic concept is to relate the Bethe-
Salpeter amplitude to the two-body propagator (four-point propagator) for
which an integral equation can be derived from perturbation theory
Feynman:1949hz ; Salpeter:1951sz ; DL . In the unitary chiral approach, the
bound state and resonance can be expressed as pole singularities in the
scattering amplitude. The bound state and resonance poles appear in the first
and second Riemann sheets, respectively Hyodo:2011ur ; Oset:1997it . In the
Bethe-Salpeter technique, the bound state gives rise to a pole in the Fourier
transform of the two-body propagator and we cannot get the continuum (or
scattering) state by solving the homogeneous Bethe-Salpeter equation DL . Both
of the two approaches include a loop integration in momentum space, which
might lead to an ultraviolet divergence. The unitary chiral approach deals
with the divergence in the so-called on-shell factorization approach
Hyodo:2011ur ; Oset:1997it ; Oller:1997ti . In the Bethe-Salpeter technique,
this divergence can be avoided which will be shown below Guo:1996jj ;
Guo:2007mm ; Xie:2010zza ; Feng:2011zzb . Besides, we solve Bethe-Salpeter
amplitude from the Bethe-Salpeter equation, while the scattering amplitude is
solved in the unitary chiral approach. The dependence of the Bethe-Salpeter
amplitude on the momentum transform reveals the structure of the bound states
Guo:1996jj ; Guo:2007mm ; Xie:2010zza ; Feng:2011zzb .

The Bethe-Salpeter techniques were developed by Feynman, Salpeter, and Bethe
Feynman:1949hz ; Salpeter:1951sz ; DL . It has been applied to theoretical
studies concerning heavy baryons and molecular bound states Guo:1996jj ;
Guo:2007mm ; Xie:2010zza ; Feng:2011zzb . In previous studies, the possible
bound states of $K\bar{K}$, $DK$, and $B\bar{K}$ have been investigated in the
Bethe-Salpeter formalism with the kernel introduced by ChPT in the ladder and
instantaneous approximations Guo:2007mm ; Xie:2010zza ; Feng:2011zzb . We will
try to study the $K^{-}p$ bound state in this framework in the present paper.
We will take these two approximations into account and consider the
interaction kernel provided by the leading order of ChPT. We will investigate
whether the bound state exists or not and study its decay in this picture. We
will also discuss the extent to which the $K^{-}p$ component contributes to
the observed $\Lambda(1405)$ resonance.

The remainder of this paper is organized as follows. In Sec. II, we derive the
Bethe-Salpeter equation for the $K^{-}p$ system in detail and present the
normalization condition of the Bethe-Salpeter amplitude. In Sec. III, the
decay of the $K^{-}p$ bound state to $\Sigma^{+}\pi^{-}$ is discussed. The
numerical results are presented in Sec. IV. In the last section, we give a
summary and some discussions.

##  II Bethe-Salpeter formalism for the bound state containing the proton and
antikaon

###  II.1 Bethe-Salpeter equation for the $K^{-}p$ system

In this section, we will derive the Bethe-Salpeter equation for the $K^{-}p$
system. We assume that the bound state exists and its mass is $M$. We denote
it by $\Lambda^{*}$. In this picture, the Bethe-Salpeter amplitude can be
defined as DL ; Guo:1996jj ; Guo:2007mm ; Xie:2010zza ; Feng:2011zzb ;
Zhang:2013gqa ; Liu:2015qfa ; Liu:2016wzh :

$\displaystyle\chi(x_{1},x_{2},P)$ $\displaystyle=$ $\displaystyle\langle
0|T\psi(x_{1})\phi(x_{2})|\Lambda^{*}\rangle,$ (1)

with $\psi(x_{1})$ and $\phi(x_{2})$ being field operators of the proton and
$K^{-}$, respectively, and $P$ being the momentum of the system. In the
momentum space, the Bethe-Salpeter amplitude, $\chi_{P}(p)$, is related to
$\chi(x_{1},x_{2},P)$ through the following equation DL :

$\displaystyle\chi(x_{1},x_{2},P)=\mathrm{e}^{\mathrm{i}PX}\int\frac{\mathrm{d}^{4}p}{(2\pi)^{4}}\chi_{P}(p)\mathrm{e}^{\mathrm{i}px},$
(2)

where $p$ and $x(=x_{1}-x_{2})$ are the relative momentum and the relative
coordinate of two constituents, respectively, and $X$ is the center of mass
coordinate which is defined as $X=\lambda_{1}x_{1}+\lambda_{2}x_{2}$, where
$\lambda_{1}=\frac{m_{1}}{m_{1}+m_{2}}$,
$\lambda_{2}=\frac{m_{2}}{m_{1}+m_{2}}$, with $m_{1}$ and $m_{2}$ being the
masses of the proton and the $K^{-}$ meson, respectively. The momentum of the
proton is $p_{1}=\lambda_{1}P+p$ and that of $K^{-}$ is
$p_{2}=\lambda_{2}P-p$. The derivation of the Bethe-Salpeter formalism for the
two fermion systems can be found in the textbook DL . In the same way, one can
prove that the form of the Bethe-Salpeter equation is still valid for the
fermion and scalar object system Guo:1996jj . Therefore, the Bethe-Salpeter
amplitude in our case satisfies the follow homogeneous integral equation DL ;
Guo:1996jj ; Guo:2007mm ; Xie:2010zza ; Feng:2011zzb ; Zhang:2013gqa ;
Liu:2015qfa ; Liu:2016wzh :

$\displaystyle\chi_{P}(p)$ $\displaystyle=$ $\displaystyle
s_{F}(\lambda_{1}P+p)\int\frac{\mathrm{d}^{4}q}{(2\pi)^{4}}K(P,p,q)\chi_{P}(q)s_{B}(\lambda_{2}P-p),$
(3)

where $s_{F}$ and $s_{B}$ are propagators of the proton and $K^{-}$,
respectively. We also define the relative longitudinal momentum $p_{l}(=v\cdot
p)$ and transverse momentum $p_{t}[=p-(v\cdot p)v]$ with $v(=P/M)$ being the
4-velocity of the bound state, and $K(P,p,q)$ is the interaction kernel that
can be described by the sum of all the irreducible graphs which cannot be
split as into two pieces by cutting two-particle lines as defined in Ref. DL .
For the propagators, we have

$\displaystyle
s_{F}(\lambda_{1}P+p)=\frac{\mathrm{i}\left[(\lambda_{1}M+p_{l})v\\!\\!\\!/+p\\!\\!\\!/_{t}+m_{1}\right]}{(\lambda_{1}M+p_{l}-\omega_{1}+\mathrm{i}\epsilon)(\lambda_{1}M+p_{l}+\omega_{1}-\mathrm{i}\epsilon)},$
$\displaystyle
s_{B}(\lambda_{2}P-p)=\frac{\mathrm{i}}{(\lambda_{2}M-p_{l}-\omega_{2}+\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}+\omega_{2}-\mathrm{i}\epsilon)},$
(4)

with $\omega_{1(2)}=\sqrt{m_{1(2)}^{2}-p_{t}^{2}}$.

In general, considering $v\\!\\!\\!/u(v,s)=u(v,s)$ , $\chi_{P}(p)$ can be
written as Zhang:2013gqa ; Liu:2015qfa ; Liu:2016wzh

$\displaystyle\chi_{P}(p)=(g_{1}+g_{2}\gamma_{5}+g_{3}\gamma_{5}p\\!\\!\\!/_{t}+g_{4}p\\!\\!\\!/_{t}+g_{5}\sigma_{\mu\nu}\varepsilon^{\mu\nu\alpha\beta}p_{t\alpha}v_{\beta})u(v,s),$
(5)

where $u(v,s)$ is the spinor of the bound state with helicity $s$ and
$g_{i}$($i=1,2\cdots 5$) are Lorentz-scalar functions. With the constraints
imposed by parity and Lorentz transformations, it is easy to prove that
$\chi_{P}(p)$ can be simplified as Zhang:2013gqa ; Liu:2015qfa ; Liu:2016wzh

$\displaystyle\chi_{P}(p)=[f_{1}(p)+f_{2}(p)p\\!\\!\\!/_{t}]u(v,s),$ (6)

in which $f_{1}(p)$ and $f_{2}(p)$ are two independent Lorentz-scalar
functions of $p$.

Figure 1: Feynman diagram for the Weinberg-Tomozawa interaction in ChPT.

The Bethe-Salpeter equation can be treated in the so-called ladder
approximation DL . In this approximation, $K(P,p,q)$ is replaced by its
lowest-order value DL . Here we adapt ChPT to describe the s-wave meson-baryon
interaction. For the $s$-wave amplitude, the most important piece is the
Weinberg-Tomozawa contact interaction at the lowest order $\mathcal{O}(p)$ of
ChPT Hyodo:2011ur ; Weingerg:1996 . According to ChPT, the Lagrangian for the
Weinberg-Tomozawa contact interaction as shown in Fig. 1 is Hyodo:2011ur ;
Ecker:1994gg ; Oset:1997it

$\displaystyle\mathcal{L}=\bar{\psi}\mathrm{i}\gamma^{\mu}\frac{C_{\mathrm{if}}}{4f^{2}}(\phi\partial_{\mu}\phi-\partial_{\mu}\phi\phi)\psi,$
(7)

where $C_{\mathrm{if}}$ is the isospin coefficient (i and f represent the
initial and final states, respectively) and $f$ corresponds to the meson decay
constant in the chiral limit at the tree level. Using this interaction term,
we obtain the formalism of $K(P,p,q)$ in terms of $p_{l}$ and $p_{t}$:

$\displaystyle K(P,p,q)$ $\displaystyle=$
$\displaystyle-\mathrm{i}C_{\mathrm{if}}\frac{1}{4f^{2}}(p\\!\\!\\!/_{2}+q\\!\\!\\!/_{2})$
(8) $\displaystyle=$
$\displaystyle-\mathrm{i}C_{\mathrm{if}}\frac{1}{4f^{2}}\left[2(\lambda_{2}M-p_{l})v\\!\\!\\!/-(p\\!\\!\\!/_{t}+q\\!\\!\\!/_{t})\right],$

where we use the covariant instantaneous approximation, $p_{l}=q_{l}$, in the
last line .

Now we substitute Eqs. (II.1), (6) and (8) into the Bethe-Salpeter equation
(3) and obtain the coupled integral equations about $f_{1}(p_{l},\,{p}_{t})$
and $f_{2}(p_{l},\,{p}_{t})$:

$\displaystyle
4f_{1}(p_{l},\,{p}_{t})=\frac{\mathrm{i}C_{\mathrm{if}}}{f^{2}}\int\frac{\mathrm{d}^{4}q}{(2\pi)^{4}}$
$\displaystyle\cdot\Bigg{[}\frac{\left[2(\lambda_{1}M+p_{l})(\lambda_{2}M-p_{l})-p_{t}^{2}-p_{t}\cdot
q_{t}\right]\cdot
f_{1}(q_{l},\,{q}_{t})}{(\lambda_{1}M+p_{l}-\omega_{1}+\mathrm{i}\epsilon)(\lambda_{1}M+p_{l}+\omega_{1}-\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}-\omega_{2}+\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}+\omega_{2}-\mathrm{i}\epsilon)}$
$\displaystyle\qquad-\frac{m_{1}(p_{t}\cdot q_{t}+q_{t}^{2})\cdot
f_{2}(q_{l},\,{q}_{t})}{(\lambda_{1}M+p_{l}-\omega_{1}+\mathrm{i}\epsilon)(\lambda_{1}M+p_{l}+\omega_{1}-\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}-\omega_{2}+\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}+\omega_{2}-\mathrm{i}\epsilon)}\Bigg{]},$
$\displaystyle
4p_{t}^{2}f_{2}(p_{l},\,{p}_{t})=\frac{\mathrm{i}C_{\mathrm{if}}}{f^{2}}\int\frac{\mathrm{d}^{4}q}{(2\pi)^{4}}$
$\displaystyle\cdot\Bigg{[}\frac{-m_{1}\left(p_{t}^{2}+p_{t}\cdot
q_{t}\right)\cdot
f_{1}(q_{l},\,{q}_{t})}{(\lambda_{1}M+p_{l}-\omega_{1}+\mathrm{i}\epsilon)(\lambda_{1}M+p_{l}+\omega_{1}-\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}-\omega_{2}+\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}+\omega_{2}-\mathrm{i}\epsilon)}$
$\displaystyle\quad+\frac{\left[2(\lambda_{1}M+p_{l})(\lambda_{2}M-p_{l})p_{t}\cdot
q_{t}-p_{t}^{2}p_{t}\cdot
q_{t}-p_{t}^{2}q_{t}^{2}\right]f_{2}(q_{l},\,{q}_{t})}{(\lambda_{1}M+p_{l}-\omega_{1}+\mathrm{i}\epsilon)(\lambda_{1}M+p_{l}+\omega_{1}-\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}-\omega_{2}+\mathrm{i}\epsilon)(\lambda_{2}M-p_{l}+\omega_{2}-\mathrm{i}\epsilon)}\Bigg{]}.$

We define the functions
$\tilde{f}_{1(2)}(p_{t})=\int\frac{\mathrm{d}p_{l}}{2\pi}f_{1(2)}(p_{l},\,p_{t})$.
In the $\Lambda^{*}$ rest frame, one has $p_{t}=(0,\,-\vec{p_{t}})$ and
$|p_{t}|=|\vec{p}_{t}|$. Performing the integration over $p_{l}$ on both sides
through the residue theorem, we find that $\tilde{f}_{1(2)}(\vec{p}_{t})$
satisfy the coupled integral equations as follows:

$\displaystyle\tilde{f}_{1}(\vec{p}_{t})=\frac{C_{\mathrm{if}}}{4f^{2}}\int\frac{\mathrm{d}^{3}\vec{q}_{t}}{(2\pi)^{3}}$
$\displaystyle\cdot\Bigg{[}\bigg{(}\frac{2(2\lambda_{1}M+\omega_{1})(M+\omega_{1})-|\vec{p}_{t}|^{2}-\vec{p}_{t}\cdot\vec{q}_{t}}{2\omega_{1}[(M+\omega_{1})^{2}-\omega_{2}^{2}]}+\frac{2\omega_{2}(M-\omega_{2})+|\vec{p}_{t}|^{2}+\vec{p}_{t}\cdot\vec{q}_{t}}{2\omega_{2}[(M-\omega_{2})^{2}-\omega_{1}^{2}]}\bigg{)}\cdot\tilde{f}_{1}(\vec{q}_{t})$
$\displaystyle\qquad\qquad-
m_{1}\cdot\bigg{(}\frac{|\vec{q}_{t}|^{2}+\vec{p}_{t}\cdot\vec{q}_{t}}{2\omega_{1}[(M+\omega_{1})^{2}-\omega_{2}^{2}]}-\frac{|\vec{q}_{t}|^{2}+\vec{p}_{t}\cdot\vec{q}_{t}}{2\omega_{2}[(M-\omega_{2})^{2}-\omega_{1}^{2}]}\bigg{)}\cdot\tilde{f}_{2}(\vec{q}_{t})\Bigg{]},$
(11) $\displaystyle\tilde{f}_{2}(\vec{p}_{t})$ $\displaystyle=$
$\displaystyle\frac{C_{\mathrm{if}}}{4(-|\vec{p}_{t}|^{2})f^{2}}\int\frac{\mathrm{d}^{3}\vec{q}_{t}}{(2\pi)^{3}}\cdot\Bigg{[}\bigg{(}\frac{m_{1}\left(|\vec{p}_{t}|^{2}+\vec{p}_{t}\cdot\vec{q}_{t}\right)}{2\omega_{1}[(M+\omega_{1})^{2}-\omega_{2}^{2}]}-\frac{m_{1}\left(|\vec{p}_{t}|^{2}+\vec{p}_{t}\cdot\vec{q}_{t}\right)}{2\omega_{2}[(M-\omega_{2})^{2}-\omega_{1}^{2}]}\bigg{)}\cdot\tilde{f}_{1}(\vec{q}_{t})$
(12)
$\displaystyle\qquad-\bigg{(}\frac{2(2\lambda_{1}M+\omega_{1})(M+\omega_{1})\vec{p}_{t}\cdot\vec{q}_{t}-|\vec{p}_{t}|^{2}\vec{p}_{t}\cdot\vec{q}_{t}-|\vec{p}_{t}|^{2}|\vec{q}_{t}|^{2}}{2\omega_{1}[(M+\omega_{1})^{2}-\omega_{2}^{2}]}$
$\displaystyle\qquad\qquad\qquad+\frac{2\omega_{2}(M-\omega_{2})\vec{p}_{t}\cdot\vec{q}_{t}+|\vec{p}_{t}|^{2}\vec{p}_{t}\cdot\vec{q}_{t}+|\vec{p}_{t}|^{2}|\vec{q}_{t}|^{2}}{2\omega_{2}[(M-\omega_{2})^{2}-\omega_{1}^{2}]}\bigg{)}\cdot\tilde{f}_{2}(\vec{q}_{t})\Bigg{]},$

with $\omega_{1(2)}=\sqrt{m_{1(2)}^{2}+|\vec{p}_{t}|^{2}}$. These two
equations involve the integrations of $q_{t}$. They look like divergent
integrations since $q_{t}$ varies from 0 to $+\infty$. However, the Bethe-
Salpeter amplitudes ($f_{1}$ and $f_{2}$) decrease to zero rapidly at the
large momentum transfer and thus there is no divergence in practice DL ;
Guo:1996jj ; Guo:2007mm ; Xie:2010zza ; Feng:2011zzb ; Feng:2011zzb .

###  II.2 Normalization condition of the Bethe-Salpeter amplitude

In Eqs. (II.1) and (II.1), we leave the normalization of $\tilde{f}_{1}$ and
$\tilde{f}_{2}$ undetermined. Following Ref.DL , the normalization condition
for the Bethe-Salpeter equation can be written as

$\displaystyle\frac{\mathrm{i}}{(2\pi)^{4}}\int\mathrm{d}^{4}p\mathrm{d}^{4}q\bar{\chi}_{P}(p)\frac{\partial}{\partial
P_{0}}[I(p,\,q,\,P)+K(p,\,q,\,P)]\chi_{P}(q)=2P_{0},$ (13)

where $I(p,\,q,\,P)$ is the inverse of the four-point propagator

$\displaystyle
I(p,\,q,\,P)=\delta^{(4)}(p-q)[s_{F}(\lambda_{1}P+p)]^{-1}[s_{B}(\lambda_{2}P-p)]^{-1}.$
(14)

One can recast the normalization condition for the Bethe-Salpeter amplitude
into the form Liu:2015qfa

$\displaystyle-\int\frac{\mathrm{d}^{4}p}{(2\pi)^{4}}\Big{\\{}\mathrm{Tr}[\alpha_{P}(p)\beta_{P}(p)s_{F}(p_{1})(\lambda_{1}\varepsilon\\!\\!\\!/)s_{F}(p_{1})s_{B}(p_{2})]$
$\displaystyle\qquad\qquad\quad+\mathrm{Tr}[\alpha_{P}(p)\beta_{P}(p)(2\lambda_{2}p_{2}\cdot\varepsilon)s_{F}(p_{1})s_{B}(p_{2})s_{B}(p_{2})]\Big{\\}}=2P_{0},$
(15)

where $\varepsilon=(1,\,0,\,0,\,0)$ and
$\alpha_{P}(p_{l},\,p_{t})$[$\beta_{P}(p_{l},\,p_{t})$] is the transverse
projection of the Bethe-Salpeter amplitude given by

$\displaystyle\alpha_{P}(p_{l},\,{p_{t}})$ $\displaystyle=$
$\displaystyle-\mathrm{i}s_{F}(p_{1})^{-1}\chi_{P}(p_{l},{p_{t}})s_{B}(p_{2})^{-1},$
$\displaystyle\beta_{P}(p_{l},\,{p_{t}})$ $\displaystyle=$
$\displaystyle-\mathrm{i}s_{B}(p_{2})^{-1}\bar{\chi}_{P}(p_{l},{p_{t}})s_{F}(p_{1})^{-1}.$
(16)

Substituting Eq. (3) into the above equations, we obtain

$\displaystyle\alpha_{P}(p_{l},\,p_{t})$ $\displaystyle=$
$\displaystyle[\tilde{h}_{1}(p_{t})+p\\!\\!\\!/_{t}\tilde{h}_{2}(p_{t})]u(v,s),$
$\displaystyle\beta_{P}(p_{l},\,p_{t})$ $\displaystyle=$
$\displaystyle\bar{u}(v,s)[\tilde{h}_{1}(p_{t})+p\\!\\!\\!/_{t}\tilde{h}_{2}(p_{t})],$
(17)

with

$\displaystyle\tilde{h}_{1}(p_{t})$ $\displaystyle=$
$\displaystyle\int\frac{\mathrm{d}^{3}q_{t}}{(2\pi)^{3}}\frac{C_{\mathrm{if}}(p_{t}\cdot
q_{t}+|{q}_{t}|^{2})\tilde{f}_{2}(q_{t})}{4f^{2}|{p}_{t}|^{2}},$
$\displaystyle\tilde{h}_{2}(p_{t})$ $\displaystyle=$
$\displaystyle\int\frac{\mathrm{d}^{3}q_{t}}{(2\pi)^{3}}\frac{C_{\mathrm{if}}(p_{t}\cdot
q_{t}+|{q}_{t}|^{2})\tilde{f}_{1}(q_{t})}{4f^{2}}.$ (18)

Then, we substitute Eqs. (II.2) and (II.2) into Eq. (II.2) and integrate out
the relative longitudinal momentum $q_{l}$. In the $\Lambda^{*}$ rest frame,
the normalization condition can be written in the following form

$\displaystyle-\int\frac{\mathrm{d}^{3}\vec{p}_{t}}{2(2\pi)^{3}}\frac{4M}{(M+\omega_{1}-\omega_{2})^{2}}\Bigg{\\{}\frac{\lambda_{1}}{2\omega_{2}(M-\omega_{1}-\omega_{2})^{2}}\Big{[}-4\tilde{h}_{1}(\vec{p}_{t})\tilde{h}_{2}(\vec{p}_{t})\vec{p}_{t}^{2}(\lambda_{2}M-\omega_{2})$
$\displaystyle\qquad\qquad\quad-\tilde{h}_{1}^{2}(\vec{p}_{t})\big{(}\lambda_{2}^{2}M^{2}+m_{1}^{2}-2\lambda_{2}M\omega_{2}+2\lambda_{2}Mm_{1}-2m_{1}\omega_{2}+\omega_{2}^{2}-\vec{p}_{t}^{2}\big{)}$
$\displaystyle\qquad\qquad\quad+\tilde{h}_{2}^{2}(\vec{p}_{t})\vec{p}_{t}^{2}\big{(}\lambda_{2}^{2}M^{2}-2\lambda_{2}Mm_{1}+m_{1}^{2}+2m_{1}\omega_{2}+\omega_{2}^{2}-\vec{p}_{t}^{2}\big{)}\Big{]}$
$\displaystyle\qquad+\frac{\lambda_{2}(M+\lambda_{1}\lambda_{2}M-\omega_{1})}{\omega_{1}(M+\omega_{1}+\omega_{2})^{2}}\Big{[}\tilde{h}_{1}^{2}(\vec{p}_{t})\big{(}m_{1}-\lambda_{1}M-\omega_{1}\big{)}+\tilde{h}_{1}(\vec{p}_{t})\tilde{h}_{2}(\vec{p}_{t})\vec{p}_{t}^{2}$
$\displaystyle\qquad\qquad\quad+\tilde{h}_{2}^{2}(\vec{p}_{t})\vec{p}_{t}^{2}\big{(}\lambda_{1}M+\omega_{1}+m_{1}\big{)}\Big{]}\Bigg{\\}}=2M.$
(19)

##  III The decay width of the $K^{-}p$ bound state

Figure 2: Feynman diagram for the decay $\Lambda^{*}\to\Sigma^{+}\pi^{-}$.

In this section, we will proceed to apply the Bethe-Salpeter technique to
derive the decay width of the $K^{-}p$ bound state. According to the
experiments, $\Lambda(1405)$ exclusively decays into $\Sigma\pi(I=0)$. In
fact, Dalitz and Deloff analyzed the $\Sigma^{+}\pi^{-}$ spectrum to extract
the mass and width of $\Lambda(1405)$ Dalitz:1991sq . This clear spectrum is
frequently shown as a representative of the $\Lambda(1405)$ spectrum and used
for the input of theoretical models Hyodo:2011ur . Therefore, we will study
the decay width of the $K^{-}p$ bound state into the above final state. As
shown in Fig. 2, the relevant interaction vertex is given in Eq. (7) at the
lowest order $\mathcal{O}(p)$ through ChPT Hyodo:2011ur ; Oset:1997it . We
define $p_{a}[=(E_{a},\,-\vec{p}_{a})]$ and $p_{b}[=(E_{b},\,-\vec{p}_{b})]$
to be the momenta of $\Sigma^{+}$ and $\pi^{-}$, respectively.
$p^{\prime}(=\lambda^{\prime}_{2}p_{a}-\lambda^{\prime}_{1}p_{b})$ is defined
as the relative momentum of $\Sigma^{+}$ and $\pi^{-}$ where
$\lambda^{\prime}_{1}=\frac{m_{a}}{m_{a}+m_{b}}$,
$\lambda^{\prime}_{2}=\frac{m_{b}}{m_{a}+m_{b}}$, with $m_{a}$ and $m_{b}$
being the masses of $\Sigma^{+}$ and $\pi^{-}$, respectively, and
$p^{\prime}_{l}(=v\cdot p^{\prime})$ is the relative longitudinal momentum of
$\Sigma^{+}$ and $\pi^{-}$. According to the kinematics of the two-body decay,
in the rest frame of the bound state one has

$\displaystyle E_{b}$ $\displaystyle=$
$\displaystyle\frac{M^{2}-m_{a}^{2}+m_{b}^{2}}{2M},\qquad
E_{a}=\frac{M^{2}-m_{b}^{2}+m_{a}^{2}}{2M},$ $\displaystyle|\vec{p}_{a}|$
$\displaystyle=$
$\displaystyle|\vec{p}_{b}|=\frac{\sqrt{(M^{2}-(m_{a}+m_{b})^{2})(M^{2}-(m_{a}-m_{b})^{2})}}{2M}.$
(20)

The differential decay width reads

$\displaystyle\mathrm{d}\Gamma=\frac{1}{32\pi^{2}}|\mathcal{M}|^{2}\frac{|\vec{p}_{a}|}{M^{2}}\mathrm{d}\Omega,$
(21)

where $\Omega$ is the solid angle of $\Sigma^{+}$.

Figure 3: The truncated Bethe-Salpeter irreducible part
$T(x_{2},\,y_{2};\,x_{1},\,y_{1})$.

Next, we will present the amplitude based on the techniques in the textbook DL
. According to the LSZ reduction formula, the S-matrix element for this
process is

$\displaystyle\langle\Sigma(p_{a})\pi(p_{b})|\Lambda^{*}(P)\rangle_{H}$
$\displaystyle=$
$\displaystyle-\text{i}\int\text{d}^{4}x_{2}\text{d}^{4}y_{2}\mathrm{e}^{\mathrm{i}p_{a}x_{2}}\mathrm{e}^{\mathrm{i}p_{b}y_{2}}\bar{u}_{\Sigma}(p_{a})$
(22)
$\displaystyle\qquad\cdot(m_{b}^{2}-\partial_{y_{2}}^{2})(\mathrm{i}\partial\\!\\!\\!/_{x_{2}}-m_{a})\langle
0|\psi(x_{2})\phi(y_{2})|\Lambda^{*}(P)\rangle_{I},$

where $I$ and $H$ represent the interaction and Heisenberg pictures,
respectively. Following Eq. 9(77) in Ref. DL , one has

$\displaystyle\langle
0|\psi(x_{2})\phi(y_{2})|\Lambda^{*}(P)\rangle_{I}=-\int\mathrm{d}^{4}x_{1}\mathrm{d}^{4}y_{1}T(x_{2},\,y_{2};x_{1},\,y_{1})\chi_{P}(x_{1},\,y_{1}),$
(23)

in which $T(x_{2},\,y_{2};\,x_{1},\,y_{1})$ is defined to be the truncated
Bethe-Salpeter irreducible part. By ’truncation’, we mean the removal of the
two propagators corresponding to the incoming lines as shown in Fig. 3.
$T(x_{2},\,y_{2};\,x_{1},\,y_{1})$ can be evaluated perturbatively:

$\displaystyle T(x_{2},\,y_{2};\,x_{1},\,y_{1})$ (24) $\displaystyle=$
$\displaystyle\frac{C_{\mathrm{if}}}{4f^{2}}\int\mathrm{d}^{4}x_{i}s_{1}(x_{i}-x_{2})s_{2}(x_{i}-y_{2})\gamma^{\mu}(\partial_{x_{i}}-\overleftarrow{\partial}_{x_{i}})\delta^{(4)}(x_{1}-x_{i})\delta^{(4)}(y_{1}-x_{i})$
$\displaystyle=$
$\displaystyle\frac{C_{\mathrm{if}}}{4f^{2}}\int\mathrm{d}^{4}x_{i}\int\frac{\mathrm{d}^{4}p_{A}}{(2\pi)^{4}}\mathrm{e}^{\mathrm{i}p_{A}(x_{i}-x_{2})}(\mathrm{i}p\\!\\!\\!/_{A}-m_{a})^{-1}\int\frac{\mathrm{d}^{4}p_{B}}{(2\pi)^{4}}\mathrm{e}^{\mathrm{i}p_{B}(x_{i}-y_{2})}(p_{B}^{2}-m_{b}^{2})^{-1}\gamma^{\mu}(\partial_{x_{i}}-\overleftarrow{\partial}_{x_{i}})$
$\displaystyle\hskip
99.58464pt\times\delta^{(4)}(x_{1}-x_{i})\delta^{(4)}(y_{1}-x_{i}).$

Substituting Eqs. (23) and (24) into Eq. (22), we have

$\displaystyle\langle\Sigma(p_{a})\pi(p_{b})|\Lambda^{*}(P)\rangle_{H}$ (25)
$\displaystyle=$
$\displaystyle\text{i}\frac{C_{\mathrm{if}}}{4f^{2}}\int\text{d}^{4}x_{2}\text{d}^{4}y_{2}\mathrm{e}^{\mathrm{i}p_{a}x_{2}}\mathrm{e}^{\mathrm{i}p_{b}y_{2}}\bar{u}(p_{a})(\mathrm{i}\partial\\!\\!\\!/_{x_{2}}-m_{a})(m_{b}^{2}-\partial_{y_{2}}^{2})\int\mathrm{d}^{4}x_{1}\mathrm{d}^{4}y_{1}\int\mathrm{d}^{4}x_{i}\delta^{(4)}(x_{1}-x_{i})\delta^{(4)}(y_{1}-x_{i})$
$\displaystyle\qquad\cdot\int\frac{\mathrm{d}^{4}p_{A}}{(2\pi)^{4}}\mathrm{e}^{\mathrm{i}p_{A}(x_{i}-x_{2})}(\mathrm{i}p\\!\\!\\!/_{A}-m_{a})^{-1}\int\frac{\mathrm{d}^{4}p_{B}}{(2\pi^{4})}\mathrm{e}^{\mathrm{i}p_{B}(x_{i}-y_{2})}(p_{B}^{2}-m_{b}^{2})^{-1}\gamma^{\mu}(\partial_{x_{i}}-\overleftarrow{\partial}_{x_{i}})\chi_{P}(x_{1},\,y_{1})$
$\displaystyle=$
$\displaystyle\mathrm{i}\frac{C_{\mathrm{if}}}{4f^{2}}\bar{u}(p_{a})\int\mathrm{d}^{4}x_{1}\mathrm{d}^{4}y_{1}\int\mathrm{d}^{4}x_{i}\mathrm{e}^{\mathrm{i}p_{a}x_{i}}\mathrm{e}^{\mathrm{i}p_{b}x_{i}}(p\\!\\!\\!/_{2}+p\\!\\!\\!/_{b})\delta^{(4)}(x_{1}-x_{i})\delta^{(4)}(y_{1}-x_{i})\chi_{P}(x_{1},\,y_{1})$
$\displaystyle=$
$\displaystyle\mathrm{i}(2\pi)^{4}\delta(P-p_{a}-p_{b})\frac{C_{\mathrm{if}}}{4f^{2}}\int\frac{\mathrm{d}^{4}p}{(2\pi)^{4}}\bar{u}(p_{a})(p\\!\\!\\!/_{2}+p\\!\\!\\!/_{b})\chi_{P}(p),$

where we use the relation
$\chi_{P}(x_{i},x_{i})=\mathrm{e}^{\mathrm{i}Px_{i}}\int\frac{\mathrm{d}^{4}p}{(2\pi)^{4}}\chi_{P}(p)$.
According to Eqs. (6) and (25), the amplitude of the
$\Lambda^{*}\to\Sigma^{+}\pi^{-}$ process is

$\displaystyle\mathcal{M}=\mathrm{i}\frac{C_{\mathrm{if}}}{4f^{2}}\int\frac{\mathrm{d}^{3}p_{t}}{(2\pi)^{3}}\bar{u}_{\Sigma}(p_{a})[(\lambda_{2}M-p_{l}^{\prime})v\\!\\!\\!/-p\\!\\!\\!/_{t}+p\\!\\!\\!/_{b}]\cdot(\tilde{f}_{1}(p_{t})+\tilde{f}_{2}(p_{t})p\\!\\!\\!/_{t})u_{\Lambda}(P),$
(26)

in which we have considered the instantaneous approximation,
$p_{l}=p_{l}^{\prime}$, again.

In the $\Lambda^{*}$ rest frame, averaging over the spins of the initial state
and summing over the spins of the final state, the unpolarized decay width of
$\Lambda^{*}$ is

$\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle\frac{|\vec{p}_{a}|}{32\pi^{5}M^{2}}\frac{C^{2}_{\mathrm{if}}}{16f^{4}}\int\mathrm{d}|\vec{q}_{t}||\vec{q}_{t}|^{2}\int\mathrm{d}|\vec{p}_{t}||\vec{p}_{t}|^{2}$
(27)
$\displaystyle\cdot\Big{\\{}4M\tilde{f}_{1}(|\vec{p}_{t}|)\tilde{f}_{1}(|\vec{q}_{t}|)\big{[}(2E_{b}+2p_{l}^{\prime})(p_{a}\cdot
p_{b})-E_{a}m_{b}^{2}+E_{a}p_{l}^{\prime
2}+2E_{b}m_{a}p_{l}^{\prime}+m_{a}m_{b}^{2}+m_{a}p_{l}^{\prime 2}\big{]}$
$\displaystyle\qquad+4\tilde{f}_{2}(|\vec{p}_{t}|)\tilde{f}_{2}(|\vec{q}_{t}|)|\vec{p}_{t}|^{2}|\vec{q}_{t}|^{2}(E_{a}+m_{a})+4\tilde{f}_{1}(|\vec{p}_{t}|)\tilde{f}_{2}(|\vec{q}_{t}|)|\vec{q}_{t}|^{2}\big{[}(p_{a}\cdot
p_{b})+E_{a}p_{l}^{\prime}$
$\displaystyle\qquad\quad+E_{b}m_{a}+m_{a}p_{l}^{\prime}\big{]}+4\tilde{f}_{2}(\vec{p}_{t})\tilde{f}_{1}(|\vec{q}_{t}|)|\vec{p}_{t}|^{2}\big{[}(p_{a}\cdot
p_{b})+E_{a}p_{l}^{\prime}+E_{b}m_{a}+m_{a}p_{l}^{\prime}\big{]}\Big{\\}},$

with $p_{l}^{\prime}=\lambda^{\prime}_{2}E_{a}-\lambda^{\prime}_{1}E_{b}$.

##  IV Numerical result

Figure 4: Numerical result for the normalized Bethe-Salpeter amplitude of the
bound state. The units of $\tilde{f_{1}}$ and $\tilde{f_{2}}$ are 1 and
$\mathrm{MeV}^{-1}$, respectively.

In this part, we will solve the Bethe-Salpeter equation numerically and try to
search for the possible solution of the $K^{-}p$ bound state. To find out the
bound state in this system, one only needs to solve the homogeneous Bethe-
Salpeter equation. One solution corresponds to a possible bound state. Since
the Bethe-Selpeter amplitude for the ground state is in fact rotationally
invariant, $\tilde{f}_{1(2)}$ depends only on $|p_{t}|$. Generally, $|p_{t}|$
varies from 0 to $+\infty$ and $\tilde{f}_{1(2)}$ would decrease to zero when
$|p_{t}|\to+\infty$. We replace $|p_{t}|$ by the variable:

$\displaystyle|p_{t}|=\left[\epsilon+50\cdot\ln\left(1+\frac{1+t}{1-t}\right)\right]\,\mathrm{MeV},$
(28)

where $\epsilon$ is a small parameter and is introduced to avoid divergence in
numerical calculations and $t$ varies from -1 to 1. We then discretize Eqs.
(II.1) and (12) into $n$ pieces ($n$ is large enough) through the Gauss
quadrature rule. The Bethe-Salpeter amplitude can be written as $n$-dimension
vectors, $f_{1(2)}^{(n)}$. The coupled integral equations become two matrix
equations $f_{1(2)}^{n}=A_{1(2)1}^{n\times n}\cdot
f_{1}^{n}+A_{1(2)2}^{n\times n}\cdot f_{2}^{n}$ [$A$ corresponds to the
coefficients in Eqs. (II.1) and (12)]. One can obtain the numerical results of
the Bethe-Salpeter amplitude by solving the eigenvalue equation obtained from
the above two matrix equations.

In our calculation, we take the values of the parameters as $m_{1}=938$ MeV,
$m_{2}=493$ MeV K.A. . According to ChPT, the isospin coefficient
$C_{\mathrm{if}}=2$ in the $K^{-}p\to K^{-}p$ coupling process Oset:1997it ;
Hyodo:2007jq . For kaons in the meson-meson interaction, $f_{K}=1.19f_{\pi}$
($f_{\pi}=93$ MeV) and we should expect similar value here Ikeda:2012au ;
Oset:1997it ; Hyodo:2007jq . It can be seen from Eqs. (II.1) and (12) that
there is only one free parameter in our model, the mass $M$ of the possible
bound state. We vary $M$ from 1300 MeV to 1450 MeV in our calculation and find
that the nontrivial solution of the eigenvalue equation exists when $M=1422$
MeV. In other words, the proton and antikaon could form a bound state in this
region and its mass is 1422 MeV. The corresponding numerical results of the
Lorentz-scalar functions in the normalized Bethe-Salpeter amplitude,
$\tilde{f}_{1}(|p_{t}|)$ and $\tilde{f}_{2}(|p_{t}|)$, are given in Fig. 4.
One should note that the units of $\tilde{f}_{1}(|p_{t}|)$ and
$\tilde{f}_{2}(|p_{t}|)$ are 1 and $\mathrm{MeV}^{-1}$, respectively. In the
following, we will take $M$, $\tilde{f}_{1}(|p_{t}|)$ and
$\tilde{f}_{2}(|p_{t}|)$ as input when calculate the decay width of the bound
state.

Table 1: The pole position (width) of the $\Lambda(1405)$ in the chiral unitary approach including the Weinberg-Tomozawa term (WT), Weinberg-Tomozawa and Born terms (WTB), and next-to-leading-order (NLO) interaction of ChPT, respectively, within the SIDDHARTA experiment constraints. K.A.  Approach | Pole 1 (MeV) | pole 2 [MeV]
---|---|---
WT | WTB | NLO | WT | WTB | NLO
Refs. Ikeda:2011pi ; Ikeda:2012au  | 1422[16] | 1421[17] | 1424[26] | 1384[90] | 1385[105] | 1381[81]
Ref. Guo:2012vv , Fit II | - | - | 1421[19] | - | - |  1388[114]
Ref. Mai:2014xna , solution#2 | - | - | 1434[10] | - | - | 1330[56]
Ref. Mai:2014xna , solution#4 | - | - | 1429[12] | - | - | 1325[90]

Then, we apply the numerical solution of the Bethe-Salpeter amplitude to
calculate the decay width of $\Lambda^{*}\to\Sigma^{+}\pi^{-}$. We use the
following input parameters K.A. ; Oset:1997it ; Hyodo:2007jq : $m_{a}$=1189
MeV, $m_{b}=139$ MeV, and $C_{\mathrm{if}}=1$. With the parameters determined
above, the decay width of the process $\Lambda^{*}\to\Sigma^{+}\pi^{-}$ in our
calculation is 15 MeV.

Several studies point out the existence of the two-pole structure in the
region of the $\Lambda(1405)$. The main component may be the $K^{-}p$ bound
state which is narrow and stable and the other is the $\Sigma\pi$ continuum
(or scattering) state K.A. ; Ikeda:2011pi ; Ikeda:2012au ; Guo:2012vv ;
Mai:2014xna . The results of the pole structure in the unitary chiral approach
within the SIDDHARTA experiment constraints are displayed in Table 1. One can
see that our result is in agreement with the pole 1 which is just below the
$K^{-}p$ threshold. This situation supports the existence of the $K^{-}p$
bound state. Furthermore, according to the PDG, the peak and width of the
$\Lambda(1405)$ resonance is $1405.1^{+1.3}_{-1.0}\,$ and $50.5\pm 2.0$ MeV,
respectively K.A. . We can see the bound state in our calculations is located
in the range of the $\Lambda(1405)$ and its decay width is quite smaller than
that of $\Lambda(1405)$. That is to say, the $K^{-}p$ bound state does exist
and could contribute to the observed $\Lambda(1405)$, but there may be some
other structures in the observed resonance region.

##  V Summary and discussion

The Bethe-Salpeter formalism has been successfully applied in many theoretical
studies concerning heavy mesons, heavy baryons, and molecular bound states
automatically including relativistic corrections. In this paper, we studied
the possible $s$-wave molecular bound state of the $K^{-}p$ system in this
formalism. Considering the interaction kernel based on the ChPT at the leading
order, we established the Bethe-Salpeter equation for the $K^{-}p$ system in
the ladder and instantaneous approximations. Then, we discretized the integral
equations and solved the eigenvalue equation numerically. We confirmed the
existence of the $s$-wave $K^{-}p$ bound state in this formalism and obtained
its Bethe-Salpeter amplitude. We also calculated the decay width of the
$\Lambda^{*}\to\Sigma^{+}\pi^{-}$ process by using the Bethe-Salpeter
amplitude. According to our calculation, the mass of the $K^{-}p$ bound state
is compatible to that of the $\Lambda(1405)$ resonance.

In this work, we used the so-called ladder approximation. One may wonder if
this approximation is a good one since higher-order graphs could give more
important contributions than the ladder graphs. In fact, the legitimacy of the
application of the ladder approximation in the Bethe-Salpeter formalism has
been studied Guo:2007mm ; Gross:1982nz ; CI . It was shown that including only
the ladder graphs in the scalar-scalar system cannot lead to the correct one-
body limit Gross:1982nz  and gauge invariance cannot be maintained within the
ladder approximation. To solve these problems, the crossed-ladder graphs
should be included at least Gross:1982nz ; CI . However, in our case, the
interaction terms at lowest order $\mathcal{O}(p)$ of ChPT, which can leads to
crossed-ladder graphs, the Born terms, mainly contribute to the $p$-wave
interaction Weingerg:1996 . So, we can adopt the ladder approximation
legitimately in our model. Another approximation we took is the instantaneous
approximation. In this approximation, the energy exchange between the
constituents is neglected. The binding energy of the bound state can be
defined as $E_{b}=M-(m_{1}+m_{2})$. In our calculation, the binding energy is
–9 MeV. This shows that the binding of the constituent particles is weak;
hence, the exchange of energy between them can be neglected.

##  VI ACKNOWLEDGEMENTs

This work was supported by the Fundamental Research Funds for the Central
Universities (Project No. 3102017OQD052), NSFC-Yunnan United Fund (Project No.
U1302267), National Natural Science Foundation of China (Projects No.
11275025, 11575023, and 11775024), and the National Science Fund for
Distinguished Young Scholars (Project No. 31325005).

## References

  * (1) K.A. Olive et al. [PDG Collaboration], Chin. Phys. C 40, 100001 (2016).
  * (2) F.K. Guo, C. Hanhart, U.G. Meißner, Q. Wang, Q. Zhao and B.S. Zou, arXiv:1705.00141 [hep-ph].
  * (3) T. Hyodo and D. Jido, Prog. Part. Nucl. Phys. 67, 55 (2012).
  * (4) R.H. Dalitz and S.F. Tuan, Phys. Rev. Lett. 2, 425 (1959).
  * (5) R.H. Dalitz and S.F. Tuan, Annals Phys. 10, 307 (1960).
  * (6) M.H. Alston et al. , Phys. Rev. Lett. 6, 698 (1961).
  * (7) R.H. Dalitz, T.C. Wong and G. Rajasekaran, Phys. Rev. 153, 1617 (1967).
  * (8) A. Bohm and R.B. Teese, Phys. Rev. D 18, 4178 (1978).
  * (9) J.M.M. Hall, W. Kamleh, D.B. Leinweber, B.J. Menadue, B.J. Owen, A.W. Thomas and R.D. Young, Phys. Rev. Lett. 114, 132002 (2015).
  * (10) D. Jido, J.A. Oller, E. Oset, A. Ramos and U.G. Meißner, Nucl. Phys. A 725, 181 (2003).
  * (11) M. Bazzi et al. [SIDDHARTA Collaboration], Phys. Lett. B 704, 113 (2011).
  * (12) M. Bazzi et al. [SIDDHARTA Collaboration], Nucl. Phys. A 881, 88 (2012).
  * (13) Y. Ikeda, T. Hyodo and W. Weise, Phys. Lett. B 706, 63 (2011).
  * (14) Y. Ikeda, T. Hyodo and W. Weise, Nucl. Phys. A 881, 98 (2012).
  * (15) M. Mai and U.G. Meißner, Nucl. Phys. A 900, 51 (2013).
  * (16) Z.-H. Guo and J.A. Oller, Phys. Rev. C 87, 035202 (2013).
  * (17) M. Mai and U.G. Meißner, Eur. Phys. J. A 51, 30 (2015).
  * (18) J. Gasser and H. Leutwyler, Nucl. Phys. B 250, 465 (1985).
  * (19) G. Ecker, Prog. Part. Nucl. Phys. 35, 1 (1995).
  * (20) E. Oset and A. Ramos, Nucl. Phys. A 635, 99 (1998).
  * (21) R.P. Feynman, Phys. Rev. 76, 749 (1949).
  * (22) E.E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232 (1951).
  * (23) D. Lurie, Particles and Fields, Interscience Publishers, New York, 1968, Chap. 9.
  * (24) J. A. Oller and E. Oset, Nucl. Phys. A 620, 438 (1997); _Erratum:_ [Nucl. Phys. A 652, 407 (1999)].
  * (25) X.-H. Guo and T. Muta, Phys. Rev. D 54, 4629 (1996).
  * (26) X.-H. Guo and X.-H. Wu, Phys. Rev. D 76, 056004 (2007).
  * (27) Z.-X. Xie, G.-Q. Feng and X.-H. Guo, Phys. Rev. D 81, 036014 (2010).
  * (28) G.-Q. Feng, Z.-X. Xie and X.-H. Guo, Phys. Rev. D 83, 016003 (2011).
  * (29) L. Zhang and X.-H. Guo, Phys. Rev. D 87, 076013 (2013).
  * (30) Y. Liu, X.-H. Guo and C. Wang, Phys. Rev. D 91, 016006 (2015).
  * (31) L.-L. Liu, C. Wang, Y. Liu and X.-H. Guo, Phys. Rev. D 95, 054001 (2017).
  * (32) S. Weingerg, The quantum theory of the field, in: modern applications, Vol. 2, Cambridge University Press, London, 1996, Chap. 19.
  * (33) R.H. Dalitz and A. Deloff, J. Phys. G 17, 289 (1991).
  * (34) T. Hyodo and W. Weise, Phys. Rev. C 77, 035204 (2008).
  * (35) F. Gross, Phys. Rev. C 26, 2203 (1982).
  * (36) C.C. Itzykson and J.-B. Zuber, Quantum field theory, Mcgraw-Hill, New York, 1985, Vol. II, Chap. 10.