We call progressive paths and rushed paths two families of Dyck
paths studied by Asinowski and Jelínek, which have the same enumerating
sequence (OEIS entry A287709). We present a bijection proving this fact.
Rushed paths turn out to be in bijection with one-sided trees,
introduced by Durhuus and Ünel, which have an asymptotic enumeration involving
a stretched exponential. We show some asymptotic results and present a random
sampling procedure for rushed paths. We conclude by presenting several other
classes of related lattice paths and directed animals that may have similar
asymptotic properties.

§ INTRODUCTION

Dyck paths—sequences of up and down steps starting and ending at $0$
and staying at nonnegative height—are probably the most famous family of
lattice paths, one of the many combinatorial classes enumerated by the Catalan
numbers (OEIS entry A000108). The aim of this article is to study two
families of Dyck paths defined in the OEIS entry A287709 and studied by
Asinowski in relation to certain rectangulation models [1].  We
paraphrase their definitions below and illustrate them in
Figure <ref>. Throughout the paper, we call height of the
path $P$, and we denote by $h(P)$, the maximum height visited by $P$.

Let $P$ be a Dyck path of height $h$. We say that $P$ is progressive
if, for $i = 2,\dotsc,h$, it visits the height $i-1$ at least twice before the
first visit at height $i$. We say that $P$ is rushed if it starts
with $h$ up steps and then never again visits the height $h$.

Left: a progressive path of height $4$ with the first and second
visits at each height marked with a red and blue dot, respectively. Right: a
rushed path of height $4$. Note that progressive paths are allowed to visit
their maximal height multiple times, while rushed paths are not.

A rushed path can be described as a run of up steps followed by a
right-to-left culminating path, as defined by Bousquet-Mélou and Ponty
[8] (a culminating path is a path that visits only nonnegative
heights and visits its final height only once). This means that the
enumeration of rushed paths of a given height is easily derived from that of
culminating paths. The addition of a run of up steps can seem like a trivial
difference; however, when considering paths irrespective of height, it makes a
significant difference in the asymptotic enumeration.

A remarkable fact, due to Asinowski and Jelínek and stated in the OEIS entry
mentioned above, is that progressive and rushed paths, when counted according
to just length, have the same enumerating sequence (this is not true when
taking height into account). This calls for a bijection, but such a bijection
has not been found so far. We present such a bijection and explain how
it behaves with respect to the height statistic.

Rushed paths, as it turns out, have already been studied by Durhuus and Ünel
in the guise of one-sided trees, defined as plane trees with a height
equal to the length of their leftmost branch [12]; upon adding a
leaf on the leftmost branch, the contour path of such trees are rushed paths
(Figure <ref>, right). Their enumerating sequence exhibit a
remarkable asymptotic estimate of the form $4^ne^{-\nu
n^{1/3}}n^{-5/6}$. This sort of stretched exponential has attracted
recent attention in many combinatorial contexts
[16, 10, 14, 13], but it highly unusual for lattice paths,
which usually have a form $\mu^nn^\gamma$, with $\gamma$ normally $0$, $-1/2$
or $-3/2$ in one dimension [3].  Culminating paths, for their
part, asymptotically number $2^n/(4n)$ <cit.>.
Pushed Dyck paths also have the stretched exponential, but involve weights on
paths rather than being a proper subset (see [16], or
[13] for their tree avatars).

The plane trees corresponding to the two Dyck paths of
Figure <ref>. The tree on the left has the property that, for any
leaf of height $i\ge2$, there is a leaf of height $i-1$ to the left.

The paper is organized as follows. In Section <ref>, we
describe our bijection between rushed and progressive paths. In
Section <ref>, we give exact and asymptotic enumeration results.
In Section <ref>, we give a random sampling algorithm for rushed
paths. Finally, in Section <ref>, we give some perspectives,
including classes of directed animals linked to progressive paths.

§ BIJECTIONS

§.§ Culminating progressive paths

Before defining our bijection between rushed and progressive paths, we need
another bijection, which goes to Dyck paths to progressive culminating
paths (culminating paths satisfying the same constraint as progressive
paths). Let $P$ be a Dyck path of height $h$, written as a word on $\{u,d\}$.
Decompose it as:

\begin{equation} \label{f}
P = A_0u\dotsm A_{h-1}u\cdot A_{h}\cdot
dB_{h-1}\dotsm dB_0
\end{equation}

where the $A_i$'s and $B_i$'s are downward Dyck paths of height at most $i$.
Let $F(P)$ be the path:

\begin{equation} \label{F}
F(P) = A_0u\cdot dB_0uA_1u\dotsm dB_{h-1}uA_hu\text.
\end{equation}

The decompositions of $P$ and $F(P)$ are illustrated in Figure <ref>
(the factors $A_0$ and $B_0$ are not shown since they are always empty).

Left: a Dyck path of height $3$ decomposed as in (<ref>). Right:
its image by $F$, a progressive culminating path of height $4$ decomposed as
in (<ref>).

The mapping $F$ is a bijection from Dyck paths of semilength $n$ and
height $h$ to progressive culminating paths of length $2n+h+1$ and
height $h+1$.

First, we check that $F(P)$ is a culminating progressive path. By
construction, the factor $dB_{i-2}uA_{i-1}u$ goes from the first visit at
height $i-1$ to the first visit at height $i$; it never visits negative height
and visits the height $i-1$ at least twice, so $F(P)$ is culminating and

Second, the decomposition (<ref>) can be recovered from $F(P)$ by cutting at
the first and second visits at height $i$, for $i = 1,\dotsc,h-1$, which
allows us to rebuild the path $P$. This shows that $F$ is a bijection.

§.§ Progressive paths

We are now ready to define our main bijection. Let $P$ be a rushed path of
height $h$. Decompose it as:

\begin{equation} \label{g}
P = u^h\cdot dC_0\dotsm dC_{h-1}
\end{equation}

where the $C_i$'s are Dyck paths of height at most $i$. Let $m$ be the maximum
height of $C_i$ for $i = 0,\dotsc,h-1$ and let $j$ be the smallest index with
$h(C_j) = m$ (necessarily, we have $j\ge m$). Let $G(P)$ be the path:

\begin{equation} \label{G}
G(P) = F(C_j) \cdot C_0d\dotsm C_{m-1}d\cdot
uC_{m}d\dotsm uC_{j-1}d\cdot d\cdot
uC_{j+1}d\dotsm uC_{h-1}d\text.
\end{equation}

The decompositions (<ref>) and (<ref>) are illustrated in
Figure <ref>.

Above: a rushed Dyck path of height $10$ decomposed as in (<ref>),
where $m = 4$ and $j = 6$ (so $C_6$ has height exactly $4$, while $C_7$, $C_8$
and $C_9$ have height at most $4$). Below: its image by $G$, a progressive
Dyck path of height $5$ decomposed as in (<ref>).

The mapping $G$ is a bijection from Dyck paths of semilength $n$
to progressive Dyck paths of semilength $n$.

First, we check that $G(P)$ is a progressive path of height $m+1$ and
length $2n$. The factor $F(C_k)$ is a progressive culminating path of
height $m+1$ by Theorem <ref>; since the factors $C_i$ have height at
most $i$ for $i = 0,\dotsc,m-1$, at most $m-1$ for $i = m,\dotsc,k-1$ and at
most $m$ for $i = k,\dotsc,h-1$, the rest of the path never visits heights
above $m+1$. Finally, the length of $G(P)$ is $\abs{P} + \abs{F(C_j)} -
\abs{C_j} - m - 1$, which is $2n$ by Theorem <ref>.

Second, the decomposition (<ref>) and thus the path $P$ can be recovered
from the path $G(P)$ by cutting at the first visit at height $m+1$, at the
first visit at height $i$ thereafter for $i=m,\dotsc,2$, at every visit at
height $1$ thereafter before the first visit at height $0$ and at every visit
at height $0$ thereafter. This shows that $G$ is a bijection.

§ ENUMERATION AND ASYMPTOTICS

§.§ Notations and classical results

Let $p_n = r_n$ be the number of progressive or rushed paths of
semilength $n$. Let $p_{n,h}$, $r_{n,h}$ and $d_{n,h}$ be the number of
progressive paths, rushed paths and Dyck paths, respectively, of
semilength $n$ and height $h$. Also let $c_{n,h}$ be the number of culminating
paths of length $n$ and height $h$ (culminating paths have odd length when $h$
is odd).  Let $R(z) = P(z)$, $R_h(z)$, $P_h(z)$, $D_h(z)$ and $C_h(z)$ be the
generating functions of these sets, where the variable $z$ counts the
semilength (so $C_h(z)$ involves powers of $z^{1/2}$ when $h$ is odd).

Let $F_h(z)$ be the Fibonacci polynomials, defined as $F_0(z) = F_1(z) = 1$
and $F_h(z) = F_{h-1}(z) - zF_{h-2}(z)$ for $h\ge2$. Also let $q := q(z)$ be
the generating function of the Catalan numbers, satisfying $q = z(1+q)^2$.

We now state some classical enumeration results for Dyck paths of bounded
height and culminating paths. First, we have:

\begin{equation} \label{DC}
D_h(z) = \frac{z^{h}}{F_h(z)F_{h+1}(z)}
\quad\text{and}\quad
C_h(z) = \frac{z^{h/2}}{F_h(z)}
\text.
\end{equation}

We can get explicit forms of these generating functions in terms of the
generating function $q := q(z)$ of Catalan numbers satisfying $q = z^2(1+q)^2$
using the formula:

\begin{equation} \label{fibo}
F_h(z) = \frac{1 - q^{h+1}}{(1-q)(1+q)^h}
\text.
\end{equation}

We also have a closed form on the number of culminating paths of height $h$.
Shifting the index $h$ to get a nicer form, we have, if $n$ and $h-1$ have the
same parity:

\begin{equation} \label{c}
c_{n,h-1} = \frac{2^{n+1}}{h}
\sum_{j=1}^{\left\lfloor\!\frac{h-1}{2}\!\right\rfloor}
\sin^2\frac{j\pi}{h} \cos^{n-1}\frac{j\pi}{h}
\text.
\end{equation}

The formula for $D_h(z)$ can be found in [11], while the one for
$C_h(z)$ is in [8] (with $q = U^2$ in Bousquet-Mélou and Ponty's
notation). The formula (<ref>) can be derived with a simple fraction
decomposition as explained in [11]; alternatively, it is a
consequence of Laplace's formula for the gambler's ruin problem [15].

§.§ Rushed paths

As stated in Section <ref>, the results of [12] on
one-sided trees translate to rushed paths. Accordingly, the results below are
essentially specializations of these results (with $\mu = 0$ in the notation
of the paper to get the uniform measure).

The generating function of rushed paths of height $h$ is:

\begin{equation} \label{Rh}
R_h(z) = \frac{z^{h}}{F_h(z)}\text.
\end{equation}

The generating function of rushed paths is:

\begin{equation} \label{R}
R(z) = \sum_{h\ge0}
%\biggl(\frac q{1+q}\biggr)^h\frac{1-q}{1-q^{h+1}}
\frac{(1-q)q^h}{(1-q^{h+1})(1+q)^h}
\text.
\end{equation}

The number of rushed paths of semilength $n$ and height $h-1$ is:

\begin{equation} \label{rh}
r_{n,h-1} = \frac{4^{n+1}}{2^hh}
\sum_{j=1}^{\left\lfloor\!\frac{h-1}{2}\!\right\rfloor}
(-1)^{j+1} \sin^2\frac{j\pi}{h} \cos^{2n-h}\frac{j\pi}{h}\text.
\end{equation}

The number of rushed paths of semilength $n$ satisfies, as $n$ tends to

\begin{equation} \label{r}
r_n = \lambda\mspace{1mu}4^ne^{-\nu n^{1/3}}n^{-5/6}
\bigl[1 + O\bigl(n^{-1/3}\bigr)\bigr]
\end{equation}

where $\lambda = (4\pi)^{5/6}(\log2)^{1/3}/\sqrt3$ and $\nu =

The identities (<ref>), (<ref>) and (<ref>) immediately follow from
the fact that $r_{n,h} = c_{2n-h,h}$, by the definition of rushed paths. The
estimate (<ref>) is Theorem 5.1 in [12], after accounting for the
shift in size between rushed paths and one-sided trees (the essence of the
argument is to note that $r_{n,h}$ takes its maximum when
$h\sim(2\pi^2n/\log2)^{1/3}$ and use Laplace's method to estimate $r_n$).

Finally, we state some useful limit law results. Let $\mathfrak R_n$ be a
uniformly distributed rushed path of semilength $n$. If $P$ is a Dyck path, we
call arch of $P$ the factor between two consecutive visits at
height $0$ and we denote by $a(P)$ the number of arches of $P$.

We have, as $n$ tends to infinity:

\begin{equation} \label{normal}
\frac{h(\mathfrak R_n) - \mu n^{1/3}}{\sigma n^{1/6}}\tod\mathcal N(0,1)
\end{equation}

where $\mu = (2\pi^2/\log 2)^{1/3}$ and $\sigma =
(2\pi^2)^{1/6}/\bigl((\log2)^{2/3}\sqrt3\bigr)$. We also have:

\begin{equation} \label{arches}
a(\mathfrak R_n)\tod\mathcal A
\quad\text{where}\quad
\mathbf E\bigl[u^\mathcal A\bigr] = \frac u{2-u}2\big.^{\sfrac{u}{2-u}-1}\text.
\end{equation}

The limit law (<ref>) is implied by Laplace's method in the proof of
Theorem 5.1 in [12].

To show (<ref>), we first note that number of arches in the Dyck path
translates to number of children of the root vertex in the plane tree. The
information about that number is contained in the local limit of the tree
(Proposition 5.6 in [12]). That limit shows that the root vertex
has $R$ children belonging to the “spine”, where $R-1$ is distributed like
$\mathrm{Poisson}(\log2)$, and $R$ critical Bienaymé-Galton-Watson
trees grafted in between, each accounting for $\mathrm{Geom}(1/2)$ more
children (see the illustration in Figure 3 in [12]). This gives the
probability generating function (<ref>).

Incidentally, this local limit gives us an intuitive picture of the arches
of $\mathfrak R_n$: it has $R$ large arches, looking at the limit like a
Poisson point process, and a geometric number of small arches in between each.

§.§ Progressive paths

We conclude this section with our own contributions concerning progressive
paths. By the bijection $G$, $P(z) = R(z)$ and $p_n = r_n$ satisfy (<ref>)
and (<ref>).

The generating function of progressive paths of height $h$ is:

\begin{equation} \label{Ph}
P_h(z) = \frac{z^{2h-1}}{F_{h-1}(z)F_h(z)F_{h+1}(z)}\text.
\end{equation}

Moreover, we have:

\begin{equation} \label{delta}
h(\mathfrak R_n) - h\bigl(G(\mathfrak R_n)\bigr)
\tod\mathcal D\quad\text{where}\quad
\mathbf E\bigl[u^{\mathcal D}\bigr] = \frac1u\biggl(1 +
\frac{4(u-1)}{2-u}2\big.^{\sfrac{u}{2-u} - 1}\biggr)
\text.
\end{equation}

In particular, the height of a random progressive path of semilength $n$
also has the limit law (<ref>).

Progressive paths of height $h$ are the product of a progressive culminating
path of height $h$ (generating function $z^{h/2}D_{h-1}(z)$ by
Theorem <ref>) and a downward culminating path of height $h+1$, minus
the last downstep (generating function $C_{h+1}(z)/z^{1/2}$). This gives

It remains to prove that $\Delta_n := h(\mathfrak R_n) - h(G(\mathfrak R_n))$
tends to the limit law $\mathcal D$. If $h(\mathfrak R_n) = h$, decompose
$\mathfrak R_n$ as in (<ref>); according to Theorem <ref>, we have
$\Delta_n \ge k$ if and only if all the $A_i$'s have height at most $h-k-1$.
If this is the case, consider the path:

\[u^{h-k+1}\cdot dC_0\dotsm dC_{h-k-1}\cdot d\cdot uC_{h-k}d\dotsm uC_{h-1}d\text.\]

This is a rushed path of length $2n+2$ and height $h-k+1$ with $a+1$ arches,
and the transformation is bijective. Denoting by $r_{n,h,a}$ the number of
rushed paths with semilength $n$, height $h$ and $a$ arches, we have

\[\mathbf P\bigl(\Delta_n \ge k \mathbin{\big\vert} h(\mathfrak R_n) = h\bigr)
= \frac{r_{n+1,h-k+1,k+1}}{r_{n,h}}\text.\]

Fix $k$ and let $n$ tend to infinity. With (<ref>), expanding the sine and
cosine, we get:

\[r_{n,h} \sim
\frac{\pi^24^{n+1}}{2^hh^3}e^{-\sfrac{\pi^2n}{h^2}}\text.\]

Since $h\sim\mu n^{1/3}$ with high probability from (<ref>), this shows
that we have the form $r_{n,h} \sim \kappa 4^ne^{-\nu n^{1/3}}n^{-1}$
for some constant $\kappa$ and therefore that the ratio
$r_{n+1,h-k+1}/r_{n,h}$ tends to $4$ with high probability. Therefore, we have
from (<ref>):

\[\mathbf P(\Delta_n \ge k) \to
4\mspace{1mu}\mathbf P(\mathcal A = k+1)\text.\]

This shows convergence in distribution to $\mathcal D$.

§ RANDOM SAMPLING

In this section, we consider the problem of finding an algorithm returning a
uniformly distributed rushed path of semilength $n$ (a random progressive path
can then simply be obtained using the bijection $G$).
The best algorithms for sampling lattice paths are usually based on
anticipated rejection: paths are drawn randomly one step at a time and those
paths that do not satisfy the wanted constraints, like staying above the
$x$-axis, are rejected during the generation. This is the idea behind the
Florentine algorithms [4, 5] and the algorithms in
<cit.>. In favorable cases, such algorithms often
achieve linear complexity in the size of the output [2].

In our case, however, using such a scheme seems difficult. If we sample a
rushed path of semilength $n$ and height $h$ by outputting $h$ up steps
followed by a random path, the probability of sampling any given path without
rejecting is $2^{-2n+h}$, but we need a probability independent of $h$ to
ensure uniformity. Moreover, a random path of length $2n$ will have a height
$O(\sqrt n)$, but we need it to be $O(n^{1/3})$. This seems to show that the
rejection rate will be prohibitive.

We instead present an algorithm based on an unranking procedure: given
an integer $a$, we find the $a$th rushed path of semilength $n$ (for some
fixed ordering of the paths). This can be easily turned into a random sampling
algorithm by drawing a random rank first. Our algorithm is loosely based on
the recursive method [19], with a subquadratic complexity and without
the need to store large tables of numbers.

§.§ Rushed paths of fixed height

Let $\mathcal D_{n,i,j,h}$ be the set of paths of length $n$, starting at
height $i$, ending at $j$ and staying at heights from $0$ to $h-1$. By
considering the height at the middle point of the path, we have when $n\ge1$:

\begin{equation} \label{D}
\mathcal D_{n,i,j,h} = \bigcup_{k=0}^{h-1}
\mathcal D_{\lfloor\!\frac n2\!\rfloor,i,k,h}
\mathcal D_{\lfloor\!\frac{n+1}2\!\rfloor,k,j,h}\text.
\end{equation}

From this identity, we can find an unranking procedure for $\mathcal
D_{n,i,j,h}$ in a divide-and-conquer fashion. This extends to rushed paths of
semilength $n$ and height $h$, which satisfy $\mathcal R_{n,h} = u^hd\mathcal

Let $d_{n,i,j,h} = \#\mathcal D_{n,i,j,h}$. It is known
<cit.> that, if $i\le j$:

\begin{equation} \label{d}
d_{n,i,j,h} = \Bigl[z^{\sfrac{n-j+i}{2}}\Bigr]\frac{F_i(z)F_{h-1-j}(z)}{F_h(z)}
\text,
\end{equation}

while the Fibonacci polynomials are given by:

\begin{equation} \label{fiboGF}
F_h(z) = \bigl[y^h\bigr]\frac1{1 - y + yz}\text.
\end{equation}

Together with the symmetry relation $d_{n,i,j,h} = d_{n,j,i,h}$, this lets us
compute these numbers, for instance with Bostan and Mori's algorithm
[6]. We are now ready to present our unranking algorithm.

$n = 1$
$j = i+1$$u$
$j = i-1$$d$

$k=0$ $h-1$
$c_1\gets d_{\lfloor\!\frac n2\!\rfloor,i,k,h}$
$c_2\gets d_{\lfloor\!\frac{n+1}2\!\rfloor,k,j,h}$
$a < c_1c_2$
$P_1\gets (a \mathbin{\mathrm{mod}} c_1\mspace{-1mu})$-th path of $\mathcal D_{\lfloor\!\frac n2\!\rfloor,i,k,h}$
$P_2\gets \lfloor a/c_1\mspace{-2mu}\rfloor$-th path of $\mathcal D_{\lfloor\!\frac{n+1}2\!\rfloor,k,j,h}$

$a\gets a - c_1c_2$

$a$th path of $\mathcal D_{n,i,j,h}$

Algorithm <ref> outputs the $a$th path of $\mathcal R_{n,h}$ in time

Here, we use the notation $\Otilde$ to mean “$O$ up to logarithmic

The correction of the algorithm is a direct consequence of (<ref>). To prove
the complexity, we show that we can compute the number $d_{n,i,j,h}$ in time
$\Otilde(nh)$. We use the algorithm in [6], which, given two
polynomials $P$ and $Q$ of degree $d$, computes $[x^n] P(x)/Q(x)$ using
$O(\log n)$ multiplications of polynomials of degree $d$. If the Taylor
coefficients have exponential growth, with fast multiplication algorithms
(e.g., the Schönage-Strassen algorithm), this can be done in
time $\Otilde(nd)$.

We distinguish two cases: if $n\ge h$, the Fibonacci polynomial (of
degree $h/2$ with coefficients exponential in $h$) can be computed from
(<ref>) in time $\Otilde(h^2)$. The number $d_{n,i,j,h}$ can then be
computed from (<ref>) in time $\Otilde(nh)$. If $n < h$, we only need the
lower $n$ coefficients of the Fibonacci polynomials. By performing all
computations only up to degree $n$, we limit the complexity to $\tilde O(nh)$.

Thus, the worst-case complexity of Algorithm <ref> is $\Otilde(nh^2)$
outside of the recursive calls. We conclude using the Master Theorem.

As written, Algorithm <ref> is fairly simple and achieves the desired
complexity without the need for any precomputation. However, in practice, one
may want to store some data to speed up computations (without changing the
dominant order of magnitude):

* the Fibonacci polynomials $F_0,\dotsc,F_h$, which take $O(h^3)$ space;

* the numbers $d_{n,i,j,h}$ for $n$ lower than a threshold $t$, with the
paths smaller than $t$ being constructed step by step with the recursive
method, which reduces the amount of recursive calls at the cost of $O(t^2h^2)$

* the numbers $d_{n,i,j,h}$ computed during the course of the algorithm,
in case the same are needed later. This takes $O(nh\log n)$ space (up to $2h$
numbers of size $n$, $4h$ of size $n/2$, etc.).

§.§ General rushed paths

To find the $a$th rushed path of semilength $n$, we need to be able to compute
the numbers $r_{n,h}$. Instead of computing all of them up front, we use a
laziness idea based on the following bound: rushed paths of height $\ge h$
start with $h$ up steps and thus number at most $2^{2n-h}$. This shows that,
for $0\le m\le n$:

\begin{equation} \label{rplus}
r_n \le \sum_{h=1}^m r_{n,h} + 2^{2n-m-1}\text.
\end{equation}

This allows us to maintain an upper bound on $r_n$ than we can improve on
demand by computing more terms, like in Algorithm <ref>.

$r^+_n\gets2^{2n-1}$ and $(r_{n,h})_{1\le h\le n}\gets$ empty array

$a\gets$ random integer between $0$ and $r^+_n-1$
$b\gets a$
$h = 1$ $n$
$r_{n,h}$ is blank
compute and store $r_{n,h}$
$r^+_n\gets r^+_n - 2^{2n-h-1} + r_{n,h}$
$b\ge r^+_n$start over

$a < r_{n,h}$output the $a$th path of $\mathcal R_{n,h}$
$a\gets a - r_{n,h}$

Uniform random rushed path of semilength $n$

Algorithm <ref> outputs a uniformly distributed path of $\mathcal
R_n$ in average time $\Otilde(n^{5/3})$.

To show that the output is uniform, we note that, if $a < r_n$, the algorithm
outputs the $a$th rushed path of semilength $n$ (ranking paths by increasing
height). If $a\ge r_n$, the algorithm ultimately starts over (in the worst
case, if $a = r_n$, it will do so after computing all of the $r_{n,h}$'s).
This shows uniformity.

To show the complexity, we first estimate how many of the terms $r_{n,h}$ need
to be computed. According to (<ref>), the term $2^{2n-m}$ becomes lower
than $r_n$ at $m = O(n^{1/3})$. Since it then decreases exponentially, this
means that the algorithm will cease starting over after computing, on average,
$O(n^{1/3})$ terms.

Finally, using the algorithm of [6], computing $r_{n,h}$ takes time
$\tilde O(nh)$, so the first $O(n^{1/3})$ terms take time $\tilde O(n^{5/3})$.
We conclude with Theorem <ref>.

Sample rushed paths and their progressive images are shown in
Figure <ref>.

Left: three random rushed paths of semilength $1000$. Right: the
images of these paths by the bijection $G$. Colors are the same as in
Figure <ref>.

§ PERSPECTIVES

Many interesting problems remain on this topic. First, classes of Dyck paths
similar to progressive paths may exhibit similar asymptotical behavior: OEIS
entry A287776 describes paths one could name doubly
progressive—both left-to-right and right-to-left progressive. The
definition of progressive paths also naturally extends to unconstrained paths:
a path is progressive if, before its first visit at every
height $i\not\in\{0,1\}$, it visits twice either $i-1$ or $i+1$. There are
many possible variations and it would be interesting to see which can be
enumerated and if results like (<ref>) and (<ref>) hold.

Dyck paths are also linked to directed animals on the triangular lattice, via
bijections with heaps of dimers [18, 9]. Adopting the
convention that animals grow towards the right side, we say that a directed
animal is acute if, for every site at ordinate $y\ne0$, there are at
least two sites to the left with ordinate $y-1$ or $y+1$. Progressive paths
are in bijection with acute half-animals (all sites above the $x$-axis, see
Figure <ref>). Full acute animals seem harder to enumerate.

Left: a progressive path. Right: the corresponding acute half-animal:
for every site not on the bottom row, there are at least two sites to the left
in the row just below.

Finally, one may want to go beyond Dyck paths, to Motzkin paths, Schröder
paths, $m$-Dyck paths, or more general models [3, 7], as well as
directed animals on different lattices. It is not clear, however, which of the
many possibilities is the “right” extension of the definition of progressive
paths. The adaptation of the asymptotic results may also prove challenging.

[1]
A. Asinowski.
Private communication, 2024.

[2]
A. Bacher and A. Sportiello.
Complexity of anticipated rejection algorithms and the
Darling-Mandelbrot distribution.
Algorithmica, 75(4):812–831, 2016.

[3]
C. Banderier and P. Flajolet.
Basic analytic combinatorics of directed lattice paths.
Theoret. Comput. Sci., 281(1-2):37–80, 2002.
Selected papers in honour of Maurice Nivat.

[4]
E. Barcucci, R. Pinzani, and R. Sprugnoli.
The random generation of directed animals.
Theoret. Comput. Sci., 127(2):333–350, 1994.

[5]
E. Barcucci, R. Pinzani, and R. Sprugnoli.
The random generation of underdiagonal walks.
Discrete Math., 139(1-3):3–18, 1995.
Formal power series and algebraic combinatorics (Montreal, PQ, 1992).

[6]
A. Bostan and R. Mori.
A simple and fast algorithm for computing the $n$-th term of a
linearly recurrent sequence.
In 2021 Symposium on Simplicity in Algorithms (SOSA), pages
118–132, 2021.

[7]
M. Bousquet-Mélou.
Discrete excursions.
Sém. Lothar. Combin., 57:Art. B57d, 23, 2006/08.

[8]
M. Bousquet-Mélou and Y. Ponty.
Culminating paths.
Discrete Math. Theor. Comput. Sci., 10(2):125–152, 2008.

[9]
M. Bousquet-Mélou and A. Rechnitzer.
Lattice animals and heaps of dimers.
Discrete Mathematics, 258:235–274, 2002.

[10]
Y.-S. Chang, M. Fuchs, H. Liu, M. Wallner, and G.-R. Yu.
Enumeration of $d$-combining tree-child networks.
In AofA 2022, 2022.

[11]
N. de Bruijn, D. Knuth, and S. Rice.
The average height of planted plane trees.
In Graph Theory and Computing, pages 15–22. Academic Press,

[12]
B. Durhuus and M. Ünel.
Local limits of one-sided trees.
La Matematica, 2023.

[13]
B. Durhuus and M. Ünel.
Trees with exponential height dependent weight.
Probability Theory and Related Fields, 186:1–45, 2023.

[14]
A. Elvey Price, W. Fang, and M. Wallner.
Compacted binary trees admit a stretched exponential.
Journal of Combinatorial Theory, Series A, 177:105306, 2021.

[15]
W. Feller.
An Introduction to Probability Theory and Its Applications,
volume 1.
Wiley, 1968.

[16]
A. J. Guttmann.
Analysis of series expansions for non-algebraic singularities.
Journal of Physics A: Mathematical and Theoretical,
48(4):045209, 2015.

[17]
C. Krattenthaler.
Handbook of enumerative combinatorics, chapter Lattice path
Chapman & Hall/CRC, 2015.

[18]
G. X. Viennot.
Heaps of pieces. I. Basic definitions and combinatorial lemmas.
In Combinatoire énumérative (Montreal, Que.,
1985/Quebec, Que., 1985), volume 1234 of Lecture Notes in Math.,
pages 321–350. Springer, Berlin, 1986.

[19]
H. Wilf.
A unified setting for sequencing, ranking, and selection algorithms
for combinatorial objects.
Advances in Mathematics - ADVAN MATH, 24:281–291, 1977.


# Short cycles of random permutations with cycle weights:
point processes approach

Oleksii Galganov<EMAIL_ADDRESS>Andrii Ilienko
<EMAIL_ADDRESS>
###### Abstract

We study the asymptotic behavior of short cycles of random permutations with
cycle weights. More specifically, on a specially constructed metric space
whose elements encode all possible cycles, we consider a point process
containing all information on cycles of a given random permutation on
$\\{1,\ldots,n\\}$. The main result of the paper is the distributional
convergence with respect to the vague topology of the above processes towards
a Poisson point process as $n\to\infty$ for a wide range of cycle weights. As
an application, we give several limit theorems for various statistics of
cycles.

###### keywords:

random permutation , cycle structure , point process , Poisson convergence

###### MSC:

[2020] 60C05 , 60G55

[1]organization=Igor Sikorsky Kyiv Polytechnic Institute, addressline=Prospect
Beresteiskyi 37, postcode=03056, city=Kyiv, country=Ukraine
[2]organization=Institute of Mathematical Statistics and Actuarial Science,
University of Bern, addressline=Alpeneggstrasse 22, city=Bern,
postcode=CH-3012, country=Switzerland

##  1 Introduction

Random permutations are a classical object of combinatorial probability.
Uniform permutations (that is, uniformly distributed random elements of the
symmetric group $\SS_{n}$) have been studied since de Montmort’s matching
problem. In recent years, there has been an extensive literature on non-
uniform permutations $\sigma_{n}$ with cycle weights, which are defined by the
probability distribution

$\mathbb{P}\left\\{\sigma_{n}=\pi\right\\}=\frac{1}{h_{n}n!}\prod_{k=1}^{\infty}\theta_{k}^{C_{k}(\pi)},\qquad\pi\in\SS_{n}.$
(1)

Here $\theta_{k}$ are non-negative parameters, $C_{k}(\pi)$ stands for the
number of $k$-cycles in $\pi$, and $h_{n}$ is a normalization ensuring that
$\sum_{\pi\in\SS_{n}}\mathbb{P}\left\\{\sigma_{n}=\pi\right\\}=1$. These
permutations were introduced in Betz and Ueltschi (2011), motivated by the
theory of Bose-Einstein condensate. For other applications and connections,
see Ercolani and Ueltschi (2014) and references therein. Note that special
cases of (1) are the uniform random permutation (with $\theta_{k}=1$ for all
$k$ and $h_{n}=1$) and Ewens random permutations (with $\theta_{k}=\theta$ for
all $k$, $h_{n}=\theta^{(n)}/n!$, and $\theta^{(n)}$ standing for the rising
factorial). The latter are based on the Ewens sampling formula which was
introduced in population genetics and subsequently found numerous
applications, see Crane (2016). A rich theory of Ewens permutations was
developed in Arratia et al. (2003).

An important subject of study of random permutations is their cycle structure
and, in particular, the asymptotic statistics of short cycles (that is, cycles
of bounded length) as $n\to\infty$. It is well known that, for Ewens
permutations $\sigma_{n}$ on $\SS_{n}$, $n\geq 1$,

$\left(C_{k}(\sigma_{n}),k\geq 1\right)\xrightarrow{d}(Z_{k},k\geq 1)$ (2)

in $\mathbb{Z}_{+}^{\infty}$, where $Z_{k}$ are independent and Poisson
distributed with means $\theta/k$, see Theorem 5.1 in Arratia et al. (2003).
In the case of more general permutations with cycle weights, Corollary 2.2 in
Ercolani and Ueltschi (2014) shows that (2) remains true with independent
Poisson distributed $Z_{k}$ with means $\theta_{k}/k$ provided the stability
condition

$\lim_{n\to\infty}\frac{h_{n-1}}{h_{n}}=1$ (3)

holds. In the same paper, it is shown that (3) is satisfied for a wide range
of asymptotics of $\theta_{k}$, from sub-exponential decay to sub-exponential
growth.

What can be said about the limiting composition of short cycles themselves?
Say, for fixed points (that is, $1$-cycles), the invariance of (1) under
relabeling suggests that, conditionally on $C_{1}(\sigma_{n})=c_{1}$, the set
of fixed points of $\sigma_{n}$ is distributed as a random equiprobable
$c_{1}$-sample from $[n]\vcentcolon=\\{1,\ldots,n\\}$ without replacement.
Similar reasoning can be given for cycles of any fixed length. However, a
rigorous description of the limiting behavior is possible only within the
framework of convergence of random point measures. Using this approach, in
Section 2 we state and prove a multivariate point processes version of (2).
Its main advantage is that it allows us to easily prove further limit theorems
for different statistics of short cycles, bypassing involved combinatorial
calculations and complicated asymptotic analysis. Various examples of such
results are given in Section 3.

##  2 Preliminaries and main result

We first introduce a metric space appropriate for describing the limiting
composition of cycles. For $k\geq 1$, let

$\mathbb{X}_{k}=\left\\{\boldsymbol{x}=(x_{1},\ldots,x_{k})\in[0,1]^{k}\colon\min\\{x_{1},\ldots,x_{k}\\}=x_{1}\right\\}$
(4)

and denote by $\rho_{k}$ the Euclidean metric on $\mathbb{X}_{k}$. The last
equality in (4) is due to the fact that any element of a cycle can be regarded
as its “beginning”.

Consider now a multi-level space
$\mathbb{X}=\bigcup_{k=1}^{\infty}\mathbb{X}_{k}$ with metric given by

$\rho(\boldsymbol{x}_{1},\boldsymbol{x}_{2})=\begin{cases}\rho_{k}(\boldsymbol{x}_{1},\boldsymbol{x}_{2}),&\boldsymbol{x}_{1},\boldsymbol{x}_{2}\in\mathbb{X}_{k},\\\
\sqrt{\max\\{k_{1},k_{2}\\}},&\boldsymbol{x}_{1}\in\mathbb{X}_{k_{1}},\,\boldsymbol{x}_{2}\in\mathbb{X}_{k_{2}},\,k_{1}\neq
k_{2}.\end{cases}$

The triangle inequality holds since
$\sup_{\boldsymbol{x}_{1},\boldsymbol{x}_{2}\in\mathbb{X}_{k}}\rho_{k}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})=\sqrt{k}$.
The metric $\rho$ makes $\mathbb{X}$ a Polish space. Moreover, $\mathbb{X}$
equipped with the Borel $\sigma$-algebra $\mathcal{B}(\mathbb{X})$ turns into
a measurable space, and $B\in\mathcal{B}(\mathbb{X})$ if and only if
$B\cap\mathbb{X}_{k}\in\mathcal{B}(\mathbb{X}_{k})$ for all $k$. This makes it
possible to define a measure on
$\left(\mathbb{X},\mathcal{B}(\mathbb{X})\right)$ by

$\lambda(B)=\sum_{k=1}^{\infty}\theta_{k}\hskip
1.0pt\lambda_{k}(B\cap\mathbb{X}_{k}),\qquad B\in\mathcal{B}(\mathbb{X}),$ (5)

where $\theta_{k}$ are defined in (1) and $\lambda_{k}$ stands for the
$k$-dimensional Lebesgue measure.

Let $\delta_{\boldsymbol{x}}=\mathds{1}\\{\boldsymbol{x}\in\cdot\\}$ be the
Dirac measure at $\boldsymbol{x}$. We will focus on the limiting behavior of
random point measures on $\left(\mathbb{X},\mathcal{B}(\mathbb{X})\right)$
given by

$\Psi_{n}=\sum_{k=1}^{\infty}\;\;\sideset{}{{}^{\neq}}{\sum}_{i_{1},\ldots,i_{k}\in[n]}\delta_{\left(\frac{i_{1}}{n},\ldots,\frac{i_{k}}{n}\right)}\mathds{1}\\!\left\\{\sigma_{n}(i_{1})=i_{2},\ldots,\sigma_{n}(i_{k})=i_{1}\right\\},$
(6)

where $\sum^{\neq}$ means that the sum is taken over $k$-tuples with distinct
entries. Moreover, since $\Psi_{n}$ is considered as a measure on
$\left(\mathbb{X},\mathcal{B}(\mathbb{X})\right)$, the latter sum includes
only those tuples in which the minimum element comes first. $\Psi_{n}$ carries
all the information about the cycle structure of $\sigma_{n}$.

Recall that the vague topology on the space of locally finite measures is
generated by the integration maps
$\nu\mapsto\int_{\mathbb{X}}f\,\mathrm{d}\nu$ for all continuous functions $f$
with bounded support; see, e.g., Resnick (1987, Section 3.4) or Kallenberg
(2017, Chapter 4) for a general exposition. Denote by $\xrightarrow{vd}$ the
distributional convergence of random point measures with underlying vague
topology.

Let $\Psi$ denote the Poisson random measure on
$\left(\mathbb{X},\mathcal{B}(\mathbb{X})\right)$ with intensity measure
$\lambda$ given by (5). The following theorem can be regarded as a point
processes extension of (2).

######  Theorem 1.

Let the stability condition (3) hold. Then $\Psi_{n}\xrightarrow{vd}\Psi$ as
$n\to\infty$.

######  Remark 1.

Let $\Psi_{n}^{(k)}$ and $\Psi^{(k)}$ be the restrictions to $\mathbb{X}_{k}$
of $\Psi_{n}$ and $\Psi$, respectively. By the restriction property of Poisson
processes (see, e.g., Theorem 5.2 in Last and Penrose (2018)), $\Psi^{(k)}$
are independent homogeneous Poisson processes with intensities $\theta_{k}$.
Due to the vague continuity of the restriction mapping
$\mu\mapsto\mu\restriction_{\mathbb{X}_{k}}$, we also have
$\Psi_{n}^{(k)}\xrightarrow{vd}\Psi^{(k)}$. In particular, as
$C_{k}(\sigma_{n})=\Psi_{n}^{(k)}(\mathbb{X}_{k})$, (2) directly follows from
Theorem 1.

For the proof, we will first need an asymptotics for probabilities of cycles.

######  Lemma 2.

Let $r\geq 1$ and
$\boldsymbol{i}^{(j)}=\bigl{(}i_{1}^{(j)},\ldots,i_{k_{j}}^{(j)}\bigr{)}$,
$j\in[r]$, be disjoint integer tuples with distinct entries such that
$\boldsymbol{i}^{(j)}/n\in\mathbb{X}_{k_{j}}$. Then, under (3), we have

$\mathbb{P}\left\\{\sigma_{n}\textit{ contains the cycles
}\boldsymbol{i}^{(1)},\ldots,\boldsymbol{i}^{(r)}\right\\}\sim\frac{\prod_{j=1}^{r}\theta_{k_{j}}}{n^{k_{1}+\ldots+k_{r}}},\qquad
n\to\infty.$ (7)

######  Proof 1.

Let $\mathcal{I}$ be the set of all entries of $\boldsymbol{i}^{(j)}$,
$j\in[r]$, and $s=\\#\mathcal{I}=k_{1}+\ldots+k_{r}$, where $\\#$ stands for
the cardinality of a set. By (1), the probability in (7) equals

$\displaystyle\sum_{\tilde{\pi}\in\SS_{[n]\setminus\mathcal{I}}}\\!\mathbb{P}\left\\{\sigma_{n}=\boldsymbol{i}^{(1)}\circ\ldots\circ\boldsymbol{i}^{(r)}\circ\tilde{\pi}\right\\}$
$\displaystyle=\sum_{\tilde{\pi}\in\SS_{[n]\setminus\mathcal{I}}}\\!\frac{1}{h_{n}n!}\prod_{k=1}^{\infty}\theta_{k}^{\\#\\{j\colon
k_{j}=k\\}+C_{k}(\tilde{\pi})}$ (8)
$\displaystyle=\frac{1}{h_{n}n!}\prod_{j=1}^{r}\theta_{k_{j}}\cdot\\!\sum_{\tilde{\pi}\in\SS_{[n]\setminus\mathcal{I}}}\prod_{k=1}^{\infty}\theta_{k}^{C_{k}(\tilde{\pi})}$
$\displaystyle=\frac{h_{n-s}\hskip
1.0pt(n-s)!}{h_{n}n!}\prod_{j=1}^{r}\theta_{k_{j}},$

where the last equality follows again from (1). Hence, the claim follows from
(3).∎

######  Proof of Theorem 1 1.

Let $\langle$ mean either $($ or $[$, and the same applies to $\rangle$. By
Theorem 4.18 in Kallenberg (2017), it suffices to prove that

  1. 1.

$\lim_{n\to\infty}\mathbb{E}\hskip 1.0pt\Psi_{n}(B)=\mathbb{E}\hskip
1.0pt\Psi(B)$ for any box $B=\bigtimes_{i=1}^{k}\langle
a_{j},b_{j}\rangle\in\mathbb{X}_{k}$, $k\geq 1$,

  2. 2.

$\lim_{n\to\infty}\mathbb{P}\left\\{\Psi_{n}(U)=0\right\\}=\mathbb{P}\left\\{\Psi(U)=0\right\\}$
for any finite union $U$ of boxes from possibly different levels
$\mathbb{X}_{k}$.

Let $B$ be a box with fixed $k\geq 1$ and $B_{\neq}$ the set of points in $B$
with distinct coordinates. By (6), we have

$\mathbb{E}\hskip
1.0pt\Psi_{n}(B)=\sideset{}{{}^{\neq}}{\sum}_{\left(i_{1},\ldots,i_{k}\right)/n\in
B}\mathbb{P}\left\\{\sigma_{n}(i_{1})=i_{2},\ldots,\sigma_{n}(i_{k})=i_{1}\right\\}.$

It follows from (8) and (7) with $r=1$ that all summands on the right-hand
side are equal and asymptotically equivalent to $\theta_{k}/n^{k}$. Thus, as
$n\to\infty$,

$\mathbb{E}\hskip
1.0pt\Psi_{n}(B)\sim\frac{\theta_{k}}{n^{k}}\cdot\\#\left(B_{\neq}\cap(\mathbb{Z}^{k}/n)\right)\to\theta_{k}\hskip
1.0pt\lambda_{k}(B)=\lambda(B)=\mathbb{E}\hskip 1.0pt\Psi(B),$

which proves (i).

We now proceed to (ii). Let us fix a finite union $U$ of boxes and denote
$U_{m}=U\cap\mathbb{X}_{m}$. So, $U_{m}=\varnothing$ for $m$ greater than some
$k\geq 1$ (the maximum dimension of boxes in the union), and
$U=\bigcup_{m=1}^{k}U_{m}$. Let $\boldsymbol{i}_{m}$ and
$\boldsymbol{i}_{m}^{(\cdot)}$ stand for integer $m$-tuples with distinct
entries and $\bigcirc_{j_{m}=1}^{r_{m}}\boldsymbol{i}_{m}^{(j_{m})}$ be the
composition
$\boldsymbol{i}_{m}^{(1)}\circ\ldots\circ\boldsymbol{i}_{m}^{(r_{m})}$ of
$r_{m}$ cycles defined by such tuples. For any $R\geq 1$, by Bonferroni’s
inequality,

$\displaystyle\mathbb{P}\left\\{\Psi_{n}(U)=0\right\\}=1-\mathbb{P}\hskip
1.0pt\biggl{\\{}\bigcup_{m=1}^{k}\bigcup_{\boldsymbol{i}_{m}/n\in
U_{m}}\\{\sigma_{n}\textit{ contains the cycle
}\boldsymbol{i}_{m}\\}\biggr{\\}}$
$\displaystyle\leq\sum_{r_{1},\ldots,r_{k}=0}^{2R}(-1)^{r_{1}+\ldots+r_{k}}\sideset{}{{}^{\ast}}{\sum}_{\begin{subarray}{c}\boldsymbol{i}_{1}^{(1)}/n,\ldots,\boldsymbol{i}_{1}^{(r_{1})}/n\in
U_{1},\\\ \vbox{\linespread{0.3}\kern 2.39996pt\hbox{.}\hbox{.}\hbox{.}\kern
1.20007pt}\\\
\boldsymbol{i}_{k}^{(1)}/n,\ldots,\boldsymbol{i}_{k}^{(r_{k})}/n\in
U_{k}\end{subarray}}\mathbb{P}\left\\{\sigma_{n}\textit{ contains
}\bigcirc_{j_{1}=1}^{r_{1}}\boldsymbol{i}_{1}^{(j_{1})}\circ\ldots\circ\bigcirc_{j_{k}=1}^{r_{k}}\boldsymbol{i}_{k}^{(j_{k})}\right\\},$
(9)

where $\sum^{\ast}$ means that the sum is taken over all unordered sets of
disjoint tuples, and a similar lower bound holds with $2R$ replaced by $2R-1$.

Due to (8), all the summands in $\sum^{\ast}$ are equal and the right-hand
side of (9) becomes

$\sum_{r_{1},\ldots,r_{k}=0}^{2R}(-1)^{r_{1}+\ldots+r_{k}}\cdot\frac{h_{n-\sum_{m=1}^{k}mr_{m}}\bigl{(}n-\sum_{m=1}^{k}mr_{m}\bigr{)}!}{h_{n}n!}\cdot\frac{\prod_{m=1}^{k}\theta_{m}^{r_{m}}}{\prod_{m=1}^{k}r_{m}!}\cdot
S_{n}(r_{1},\ldots,r_{k}),$ (10)

where

$S_{n}(r_{1},\ldots,r_{k})=\sum_{\begin{subarray}{c}\boldsymbol{i}_{1}^{(1)}/n,\ldots,\boldsymbol{i}_{1}^{(r_{1})}/n\in\mathbb{X}_{1},\\\
\vbox{\linespread{0.3}\kern 2.39996pt\hbox{.}\hbox{.}\hbox{.}\kern
1.20007pt}\\\
\boldsymbol{i}_{k}^{(1)}/n,\ldots,\boldsymbol{i}_{k}^{(r_{k})}/n\in\mathbb{X}_{k}\end{subarray}}\mathds{1}\\{\textit{all
tuples are disjoint}\hskip
1.0pt\\}\prod_{m=1}^{k}\mathds{1}\bigl{\\{}\boldsymbol{i}_{m}^{(1)}/n,\ldots,\boldsymbol{i}_{m}^{(r_{m})}/n\in
U_{m}\bigr{\\}},$

and division by $\prod_{m=1}^{k}r_{m}!$ is due to the fact that the sum in the
definition of $S_{n}$ is taken over ordered sets of tuples. Note that
$\frac{S_{n}(r_{1},\ldots,r_{k})}{n^{r_{1}+\ldots+kr_{k}}}$ can be viewed as
an integral sum for

$\int_{U_{1}^{r_{1}}\times\ldots\times
U_{k}^{r_{k}}}\mathds{1}\left\\{\textit{all components of all
}\boldsymbol{x}_{m}^{(j_{m})}\textit{ are
distinct}\right\\}\;\prod_{m=1}^{k}\prod_{j_{m}=1}^{r_{m}}\mathrm{d}\boldsymbol{x}_{m}^{(j_{m})}=\prod_{m=1}^{k}\left(\lambda_{m}(U_{m})\right)^{r_{m}}.$

Hence, by (3), each summand in (10) converges as $n\to\infty$ to

$(-1)^{r_{1}+\ldots+r_{k}}\prod_{m=1}^{k}\frac{\left(\theta_{m}\lambda_{m}(U_{m})\right)^{r_{m}}}{r_{m}!}.$

It now follows from (9) and a similar lower bound that

$\displaystyle\sum_{r_{1},\ldots,r_{k}=0}^{2R-1}(-1)^{r_{1}+\ldots+r_{k}}\prod_{m=1}^{k}\frac{\left(\theta_{m}\lambda_{m}(U_{m})\right)^{r_{m}}}{r_{m}!}$
$\displaystyle\leq\lim_{n\to\infty}\mathbb{P}\left\\{\Psi_{n}(U)=0\right\\}$
$\displaystyle\leq\sum_{r_{1},\ldots,r_{k}=0}^{2R}(-1)^{r_{1}+\ldots+r_{k}}\prod_{m=1}^{k}\frac{\left(\theta_{m}\lambda_{m}(U_{m})\right)^{r_{m}}}{r_{m}!}.$

Letting $R\to\infty$ finally yields

$\displaystyle\lim_{n\to\infty}\mathbb{P}\left\\{\Psi_{n}(U)=0\right\\}$
$\displaystyle=\sum_{r_{1},\ldots,r_{k}=0}^{\infty}(-1)^{r_{1}+\ldots+r_{k}}\prod_{m=1}^{k}\frac{\left(\theta_{m}\lambda_{m}(U_{m})\right)^{r_{m}}}{r_{m}!}$
$\displaystyle=\exp\hskip
1.0pt\Bigl{\\{}-\sum_{m=1}^{k}\theta_{m}\lambda_{m}(U_{m})\Bigr{\\}}=\exp\hskip
1.0pt\\{-\lambda(U)\\}=\mathbb{P}\left\\{\Psi(U)=0\right\\}.$

This concludes the proof of (ii) and, hence, that of Theorem 1.∎

##  3 Limit theorems for statistics of short cycles

Theorem 1 allows us to derive limiting distributions for various statistics of
short cycles. In what follows, we will always assume that the stability
condition (3) is satisfied.

We first give a general result which covers the case of additive statistics.

######  Proposition 3.

Let $k\geq 1$ and $f_{m}\colon\mathbb{X}_{m}\to[0,\infty)$, $m\in[k]$, be a
family of continuous functions. Denote by $\mathcal{C}_{m}(\sigma_{n})$ the
set of all $m$-cycles in $\sigma_{n}$. Then

$\sum_{m=1}^{k}\sum_{c\hskip 1.0pt\in\hskip
1.0pt\mathcal{C}_{m}(\sigma_{n})}f_{m}\left(\frac{c}{n}\right)\xrightarrow{d}S,\qquad
n\to\infty,$ (11)

where $c\in\mathcal{C}_{m}(\sigma_{n})$ is understood as an integer tuple
whose minimum element comes first, and the limiting random variable $S$ is
defined by its Laplace transform

$\mathbb{E}\exp\hskip 1.0pt\\{-tS\\}=\exp\hskip
1.0pt\Bigl{\\{}-\sum_{m=1}^{k}\theta_{m}\int_{\mathbb{X}_{m}}\left(1-\mathrm{e}^{-tf_{m}(\boldsymbol{x})}\right)\,\mathrm{d}\boldsymbol{x}\Bigr{\\}},\qquad
t\geq 0.$ (12)

######  Proof 2.

Define $f\colon\mathbb{X}\to[0,\infty)$ by
$f(\boldsymbol{x})=\sum_{m=1}^{k}f_{m}(\boldsymbol{x})\mathds{1}\\{\boldsymbol{x}\in\mathbb{X}_{m}\\}$.
The left-hand side of (11) can be written as
$\int_{\mathbb{X}}f(\boldsymbol{x})\,\Psi_{n}(\mathrm{d}\boldsymbol{x})$, and
$f$ is continuous with bounded support. By Theorem 1 and Lemma 4.12 in
Kallenberg (2017), (11) holds with
$S=\int_{\mathbb{X}}f(\boldsymbol{x})\,\Psi(\mathrm{d}\boldsymbol{x})$. The
Laplace transform of $S$ is thus of the form

$\mathbb{E}\exp\hskip 1.0pt\\{-tS\\}=\mathbb{E}\exp\hskip
1.0pt\Bigl{\\{}-\int_{\mathbb{X}}tf(\boldsymbol{x})\,\Psi(\mathrm{d}\boldsymbol{x})\Bigr{\\}},$

which coincides with the right-hand side of (12) due to the form of the
Laplace functional of a Poisson random measure, see, e.g., Proposition 3.6 in
Resnick (1987).∎

As an example, we give a limit theorem for the sum $S_{n}^{(k)}$ of elements
in all $k$-cycles of $\sigma_{n}$.

######  Proposition 4.

  1. 1.

$\frac{S_{n}^{(k)}}{n}\xrightarrow{d}S^{(k)}$ as $n\to\infty$, where $S^{(k)}$
is defined by its Laplace transform

$\mathbb{E}\exp\left\\{-tS^{(k)}\right\\}=\exp\hskip
1.0pt\biggl{\\{}\frac{\theta_{k}}{k}\biggl{(}\left(\frac{1-\mathrm{e}^{-t}}{t}\right)^{k}-1\biggr{)}\biggr{\\}},\qquad
t>0.$ (13)

  2. 2.

If $k=1$, that is, for $S_{n}^{(1)}$ being the sum of fixed points,

$\mathbb{P}\left\\{S^{(1)}\leq
x\right\\}=\mathrm{e}^{-\theta_{1}}\sum_{j=0}^{\lfloor
x\rfloor}\frac{(-1)^{j}}{j!}(\theta_{1}(x-j))^{\frac{j}{2}}I_{j}\\!\left(2\sqrt{\theta_{1}(x-j)}\right),\qquad
x\geq 0,$ (14)

where $I_{j}$ is the modified Bessel function of the first kind, see, e.g., §
10.25 (ii) in Olver et al. (2010).

######  Proof 3.

For $\boldsymbol{x}\in\mathbb{X}_{k}$, let $f_{k}(\boldsymbol{x})$ denote the
sum of all its components. By Proposition 3, (i) holds with

$\mathbb{E}\exp\left\\{-tS^{(k)}\right\\}=\exp\hskip
1.0pt\Bigl{\\{}-\theta_{k}\int_{\mathbb{X}_{k}}\left(1-\mathrm{e}^{-tf_{k}(\boldsymbol{x})}\right)\,\mathrm{d}\boldsymbol{x}\Bigr{\\}},$

where the integral on the right-hand side, due to symmetricity of $f_{k}$,
equals

$\frac{1}{k}\int_{[0,1]^{k}}\left(1-\mathrm{e}^{-tf_{k}(\boldsymbol{x})}\right)\,\mathrm{d}\boldsymbol{x}=\frac{1}{k}\Bigl{(}1-\Bigl{(}\int_{0}^{1}\mathrm{e}^{-tx}\,\mathrm{d}x\Bigr{)}^{k}\,\Bigr{)}.$

This yields (13).

To prove (14), we first note that

$\mathbb{E}\exp\left\\{-tS^{(1)}\right\\}=\mathbb{E}\int_{0}^{\infty}t\mathrm{e}^{-tx}\,\mathds{1}\\{S^{(1)}\leq
x\\}\,\mathrm{d}x=t\int_{0}^{\infty}\mathrm{e}^{-tx}\,\mathbb{P}\hskip
1.0pt\\{S^{(1)}\leq x\\}\,\mathrm{d}x,$

cf. Liu (2020). Hence, $\mathbb{P}\left\\{S^{(1)}\leq x\right\\}$ is the
inverse Laplace transform of the right-hand side in (13) for $k=1$ multiplied
by $t^{-1}$, that is, in expanded form, of the function

$G(t)=\sum_{j=0}^{\infty}\frac{(-\theta_{1})^{j}}{j!}\mathrm{e}^{-\theta_{1}}\hskip
1.0ptt^{-j-1}\mathrm{e}^{\frac{\theta_{1}}{t}}\mathrm{e}^{-jt},\qquad t>0.$
(15)

By the time shifting property and Erdélyi et al. (1954, eq. (5.5.35)), the
inverse Laplace transform of the $j$-th summand $G_{j}(t)$ in (15) equals

$F_{j}(x)=\frac{(-1)^{j}}{j!}e^{-\theta_{1}}\hskip
1.0pt(\theta_{1}(x-j))^{\frac{j}{2}}I_{j}\\!\left(2\sqrt{\theta_{1}(x-j)}\right)\mathds{1}\\{x\geq
j\\}.$

This means that, for all $t>0$ and $R\in\mathbb{N}$,

$\int_{0}^{\infty}\\!\mathrm{e}^{-tx}\sum_{j=0}^{R}F_{j}(x)\,\mathrm{d}x=\sum_{j=0}^{R}G_{j}(t).$
(16)

Since $I_{j}(x)$ decreases in $j$ and increases in $x$ (see, e.g., § 10.37 in
Olver et al. (2010)), and $I_{0}(x)\sim\frac{\mathrm{e}^{x}}{\sqrt{2\pi x}}$
as $x\to\infty$ (ibid., § 10.30 (ii)), we have

$\Bigl{|}\sum_{j=0}^{R}F_{j}(x)\Bigr{|}\leq
e^{-\theta_{1}}I_{0}\\!\left(2\sqrt{\theta_{1}x}\right)\sum_{j=0}^{\infty}\frac{1}{j!}(\theta_{1}x)^{\frac{j}{2}}\sim
e^{-\theta_{1}}\frac{\mathrm{e}^{3\sqrt{\theta_{1}x}}}{2\pi^{\frac{1}{2}}(\theta_{1}x)^{\frac{1}{4}}},\qquad
x\to\infty,$

which makes it possible to apply dominated convergence to (16). Hence, the
inverse Laplace transform of $G$ is $\sum_{j=0}^{\infty}F_{j}$, which yields
(14).∎

We now turn to some examples of non-additive statistics. For $k\geq 2$, let us
call the range of a cycle the difference between its maximum and minimum
elements and denote by $r_{n}^{(k)}$ (resp., $R_{n}^{(k)}$) the minimum
(maximum) range among all $k$-cycles in $\sigma_{n}$. If there are no
$k$-cycles, we set $r_{n}^{(k)}=n$ and $R_{n}^{(k)}=0$.

######  Proposition 5.

$\frac{r_{n}^{(k)}}{n}\xrightarrow{d}r^{(k)}$ and
$\frac{R_{n}^{(k)}}{n}\xrightarrow{d}R^{(k)}$ as $n\to\infty$, where $r^{(k)}$
(resp., $R^{(k)}$), $k\geq 2$, have CDF’s of the form

$\displaystyle\mathbb{P}\left\\{r^{(k)}\leq x\right\\}=1-\exp\hskip
1.0pt\Bigl{\\{}-\frac{\theta_{k}}{k}\left(kx^{k-1}-(k-1)x^{k}\right)\Bigr{\\}},$
(17) $\displaystyle\mathbb{P}\left\\{R^{(k)}\leq x\right\\}=\exp\hskip
1.0pt\Bigl{\\{}\frac{\theta_{k}}{k}\left(kx^{k-1}-(k-1)x^{k}-1\right)\Bigr{\\}},$
(18)

as $x\in[0,1)$ and $0$ (resp., $1$) to the left (right ).

######  Proof 4.

For $\boldsymbol{x}\in\mathbb{X}_{k}$, let $f_{k}(\boldsymbol{x})$ denote the
difference between its maximum and minimum components. It follows from the
interpretation of vague convergence in Proposition 3.13 of Resnick (1987) that
the function which maps a finite point measure $\mu$ on $\mathbb{X}_{k}$ into
$\min_{\mu\\{\boldsymbol{x}\\}\geq 1}f_{k}(\boldsymbol{x})$ is vaguely
continuous. Hence, by Remark 1 and continuous mapping theorem,
$\frac{r_{n}^{(k)}}{n}\xrightarrow{d}r^{(k)}$ holds with
$r^{(k)}=\min_{\Psi^{(k)}\\{\boldsymbol{x}\\}\geq 1}f_{k}(\boldsymbol{x})$,
and the CDF of $r^{(k)}$ is of the form

$\displaystyle\mathbb{P}\left\\{r^{(k)}\leq x\right\\}$
$\displaystyle=1-\mathbb{P}\left\\{\Psi^{(k)}\\{\boldsymbol{x}\in\mathbb{X}_{k}\colon
f_{k}(\boldsymbol{x})\leq x\\}=0\right\\}$ (19) $\displaystyle=1-\exp\hskip
1.0pt\bigl{\\{}-\theta_{k}\hskip
1.0pt\lambda_{k}\\{\boldsymbol{x}\in\mathbb{X}_{k}\colon
f_{k}(\boldsymbol{x})\leq x\\}\bigr{\\}}.$

The value of the Lebesgue measure on the right-hand side is

$\int_{0}^{1}\biggl{(}\int_{x_{1}}^{\min\\{x_{1}+x,1\\}}\mathrm{d}x_{2}\ldots\mathrm{d}x_{k}\biggr{)}\mathrm{d}x_{1}=x^{k-1}(1-x)+\frac{x^{k}}{k},$

which together with (19) yields (17). The proof of (18) is similar.∎

As a final example, we consider some statistics of fixed points. Let $m_{n}$
be the minimum fixed point ($n+1$ in case there is none), $M_{n}$ the maximum
one ($0$ in that case), $\delta_{n}$ the minimum spacing between fixed points
(the two extreme spacings of lengths $m_{n}$ and $n+1-M_{n}$ are also taken
into account), and $\Delta_{n}$ the maximum one.

######  Proposition 6.

$\left(\frac{m_{n}}{n},\frac{M_{n}}{n},\frac{\delta_{n}}{n},\frac{\Delta_{n}}{n}\right)\xrightarrow{d}(m,M,\delta,\Delta)$
as $n\to\infty$ with

$\displaystyle\mathbb{P}\left\\{m\leq x\right\\}=1-\exp\hskip
1.0pt\\{-\theta_{1}x\\},\qquad\mathbb{P}\left\\{M\leq x\right\\}=\exp\hskip
1.0pt\\{\theta_{1}(x-1)\\},\qquad x\in[0,1),$
$\displaystyle\delta\overset{d}{=}\frac{X_{\nu+1}}{(\nu+1)\sum_{i=1}^{\nu+1}X_{i}},\qquad\Delta\overset{d}{=}\frac{\sum_{i=1}^{\nu+1}\frac{X_{i}}{i}}{\sum_{i=1}^{\nu+1}X_{i}},$

where $\nu$ is Poisson distributed with parameter $\theta_{1}$, $X_{i}$ are
exponentially distributed with unit mean, and all these are independent.

######  Proof 5.

As in the proof of Proposition 5, the convergence takes place with $m$ being
the leftmost point of $\Psi^{(1)}$, $M$ the rightmost one, $\delta$ the
minimum spacing, and $\Delta$ the maximum one. All that remains is to derive
the corresponding CDF’s and distributional equalities. For $m$ and $M$, this
follows from

$\displaystyle\mathbb{P}\left\\{m\leq
x\right\\}=1-\mathbb{P}\left\\{\Psi^{(1)}[0,x]=0\right\\}=1-\exp\hskip
1.0pt\\{-\theta_{1}x\\},$ $\displaystyle\mathbb{P}\left\\{M\leq
x\right\\}=\mathbb{P}\left\\{\Psi^{(1)}(x,1]=0\right\\}=\exp\hskip
1.0pt\\{-\theta_{1}(1-x)\\}.$

We now turn to the equalities for $\delta$ and $\Delta$. Since $\Psi^{(1)}$ is
a homogeneous Poisson process with intensity $\theta_{1}$, the conditional
distribution of $\delta$ (resp., $\Delta$) given
$\Psi^{(1)}(\mathbb{X}_{1})=r$ coincides with that of the minimum (maximum)
spacing $d_{r}$ ($D_{r}$) for $r$ independent random variables, uniformly
distributed on $[0,1]$, see, e.g., Proposition 3.8 in Last and Penrose (2018).

It follows from the theory of uniform spacings (see, e.g., Holst (1980), p.
625) that

$d_{r}\overset{d}{=}\frac{\min Y_{i}}{\sum_{i=1}^{r+1}Y_{i}},\qquad
D_{r}\overset{d}{=}\frac{\max Y_{i}}{\sum_{i=1}^{r+1}Y_{i}},$ (20)

where $Y_{1},\ldots,Y_{r+1}$ are independent $\mathsf{Exp}(1)$-distributed
random variables. Let $Y_{(i)}$, $i\in[r+1]$, be their order statistics,
$Y_{(0)}=0$, and $\tau_{i}=Y_{(i)}-Y_{(i-1)}$. By Sukhatme–Rényi decomposition
(see, e.g., Theorem 4.6.1 in Arnold et al. (2008)), $\tau_{i}$ are independent
and $\mathsf{Exp}(r-i+2)$-distributed with

$\min Y_{i}=\tau_{1},\qquad\max
Y_{i}=\sum_{i=1}^{r+1}\tau_{i},\qquad\sum_{i=1}^{r+1}Y_{i}=\sum_{i=1}^{r+1}(r-i+2)\tau_{i}.$

Denoting $X_{r-i+2}=(r-i+2)\tau_{i}\sim\mathsf{Exp}(1)$, we can rewrite (20)
as

$d_{r}\overset{d}{=}\frac{X_{r+1}}{(r+1)\sum_{i=1}^{r+1}X_{i}},\qquad
D_{r}\overset{d}{=}\frac{\sum_{i=1}^{r+1}\frac{X_{i}}{i}}{\sum_{i=1}^{r+1}X_{i}}.$

Since this holds for any $r\geq 0$, the claim follows by deconditioning.

## References

  * Arnold et al. (2008) Arnold, B.C., Balakrishnan, N., Nagaraja, H.N., 2008.  A first course in order statistics.  Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA.
  * Arratia et al. (2003) Arratia, R., Barbour, A.D., Tavaré, S., 2003.  Logarithmic combinatorial structures: a probabilistic approach.  EMS Monographs in Mathematics, European Mathematical Society (EMS), Zürich.
  * Betz and Ueltschi (2011) Betz, V., Ueltschi, D., 2011\.  Spatial random permutations with small cycle weights.  Probab. Theory Related Fields 149, 191–222.
  * Crane (2016) Crane, H., 2016.  The ubiquitous Ewens sampling formula.  Statist. Sci. 31, 1–19.
  * Ercolani and Ueltschi (2014) Ercolani, N.M., Ueltschi, D., 2014\.  Cycle structure of random permutations with cycle weights.  Random Structures Algorithms 44, 109–133.
  * Erdélyi et al. (1954) Erdélyi, A., Magnus, W., Oberhettinger, F., Tricomi, F.G., 1954\.  Tables of integral transforms. Vol. I.  McGraw-Hill Book Co., Inc., New York-Toronto-London.
  * Holst (1980) Holst, L., 1980.  On the lengths of the pieces of a stick broken at random.  J. Appl. Probab. 17, 623–634.
  * Kallenberg (2017) Kallenberg, O., 2017.  Random measures, theory and applications.  Probability Theory and Stochastic Modelling, Springer, Cham.
  * Last and Penrose (2018) Last, G., Penrose, M., 2018\.  Lectures on the Poisson process.  Institute of Mathematical Statistics Textbooks, Cambridge University Press, Cambridge.
  * Liu (2020) Liu, Y., 2020.  A general treatment of alternative expectation formulae.  Statist. Probab. Lett. 166, 108863, 8.
  * Olver et al. (2010) Olver, F.W.J., Lozier, D.W., Boisvert, R.F., Clark, C.W. (Eds.), 2010\.  NIST handbook of mathematical functions.  U.S. Department of Commerce, National Institute of Standards and Technology, Washington, DC; Cambridge University Press, Cambridge.
  * Resnick (1987) Resnick, S.I., 1987.  Extreme values, regular variation, and point processes.  Applied Probability. A Series of the Applied Probability Trust, Springer-Verlag, New York.


# Thermal equilibrium of a Brownian particle with coordinate dependent
diffusion: comparison of the Boltzmann and a modified Boltzmann distribution
with experimental results

Abhinav Dhawan  A. Bhattacharyay<EMAIL_ADDRESS>Indian
Institute of Science Education and Research, Pune, India

(June 7, 2022)

###### Abstract

In this paper we compare the Boltzmann distribution with a modified Boltzmann
distribution, that results from an Itô-process considering thermal equilibrium
of a Brownian particle with coordinate dependent diffusion, in the light of an
existing experiment. The experiment was reported in 1994 by Faucheux and
Libchaber. The experiment made use of direct tracking of diffusion of Brownian
particles near a wall. Results of this experiment allows us to compare the
Boltzmann and the modified Boltzmann distribution without making use of any
adjustable parameter. A comparison of these two distributions with the
experimental results lends support to the consideration of thermodynamic
equilibrium of a Brownian particle with coordinate-dependent diffusion to be
an Itô-process.

###### pacs:

05.10.Gg, 05.20.-y,05.40.Jc,05.70.-a

##  I Introduction

Understanding position distribution of a confined Brownian particle (BP) with
coordinate dependent diffusion and in thermal equilibrium with a heat-bath of
constant temperature is important due to immense possibility of its
application [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] across disciplines. A BP near a
wall or, similarly, near other BP’s see coordinate dependence of effective-
diffusion as a rule than an exception. The underlying reason of this
coordinate (or space) dependence of diffusion is considered mostly to be
hydrodynamic. The phrase effective-diffusion, in the present context,
indicates local renormalization of the diffusivity of a BP by the average
presence of hydrodynamic modes due to the proximity to walls and interfaces.

When it is possible to define/measure the coordinate dependence (and time
independence) of diffusivity of a BP, the consideration is that the BP, at the
least locally, has thermally equilibrated with the heat-bath. This spatially
inhomogeneous Brownian motion models a slower local isotropic transport of a
much larger BP in comparison to the time and length and scales of the bath
degrees of freedom. These bath degrees of freedom could be those faster
hydrodynamic modes whose average effect is modeled by Brownian motion. This is
the standard wisdom which even underlies the existence of coordinate-dependent
diffusivity because the local method of measurement is the same using mean-
square-displacement (MSD) against time. Due to this definition/measurement of
diffusivity the transport it represents is always locally isotropic even when
the diffusivity is coordinate-dependent.

It is important to note that, only the locally isotropic part of the average
local-transport caused by hydrodynamic effects should be modeled by a
diffusion process which has space-dependent diffusivity. The locally
anisotropic part of the average transport needs the recognition of an
effective force to account for the drift it would produce in the context of
the Langevin dynamics. If one tries to model existence of such a drift due to
hydrodynamic reasons by a noise term in the Langevin dynamics, one would end
up using correlated noise and the thermal nature of the noise will get
distorted. This anisotropic part of local transport due to hydrodynamic
effects is of no concern to us in the present context and we assume here that
the process does not involve such a local drift unaccounted for by the known
interactions. It is absolutely technically possible to keep the diffusive
transport isotropic at every point in space even when the diffusivity changes
from point to point and that is the key consideration underlying the presence
of an uncorrelated thermal noise modelling Brownian motion.

Given this context, the question that may arise is - what is the spatial
distribution of such a BP over an inhomogeneous space where the source of
inhomogeneity (space-dependent diffusivity) cannot feature in the Hamiltonian
of the system? This question has not been paid much attention to until
recently. Bhattacharyay and coworkers [11, 12, 13] have shown that, for
thermodynamic equilibrium, such a coordinate-dependent diffusion process is
essentially an Itô-process and that gives a modified Boltzmann distribution
(MBD) for thermal equilibrium. The resulting MBD of an Itô-process can be
written (for example in one-dimensional space) as

$P(x)=\frac{N}{D(x)}e^{-U(x)/[D(x)\Gamma(x)]},$ (1)

where $D(x)$ is the coordinate-dependent diffusivity and $\Gamma(x)$ is the
coordinate-dependent damping coefficient. The interaction potential that
confines the particle to equilibrate with a heat-bath is $U(x)$. Obviously, in
the above expression, $P(x)$ is the probability density at coordinate $x$
where $N$ is the normalization constant. The Langevin dynamics does not in any
way ensure local validity of the Stokes-Einstein relation
$D(x)\Gamma(x)=k_{B}T$ where $k_{B}$ is the Boltzmann constant and $T$ is the
temperature. The Stokes-Einstein relation is always imposed to the Langevin
dynamics. If one does so to an Itô-process, one would get the Boltzmann factor
appear in the expression of the MBD. It has been shown that the Fick’s law as
it comes from the Kramers-Moyal expansion is the correct one to consider [11]
in any convention (Itô or Stratonovich) and the Stokes-Einstein relation
corresponding to a constant temperature everywhere in space is consistent with
an Itô process [12].

Bhattacharyay has shown that the equilibrium distribution of a generalized
Langevin dynamics (involving inertial term and velocity-distribution thereof)
under the local validity of Stokes-Einstein relation has the over-damped limit
which converges to the equilibrium distribution of an Itô-process [12]. This
is a crucial observation in this context because an over-damped dynamics
(corresponding, for example, to an Itô-process) does not contain any
information of velocity distribution and, hence, temperature. Temperature is
imposed to an over-damped dynamics by the imposition of the Stokes-Einstein
relation. One has to, therefore, understand if at all the imposition of the
Stokes-Einstein relation is consistent with the stochastic process in its
over-damped limit. These observations have been further extended to random-
walk methods developed by Maniar and Bhattacharyay in [13] where the
maintenance of a constant temperature everywhere in space by the local-
validity of the Stokes-Einstein relation is confirmed. Such a system is seen
to correspond to an Itô-process at the continuum limit of the random walk
[13].

On the contrary, theoretical imposition of the Boltzmann distribution (BD) has
been the basis of theoretically considering equilibrium in the context of
coordinate-dependent diffusion in existing literature so far [14, 15, 16, 17,
18, 19]. In this connection, to get the Boltzmann distribution for such a BP,
one treats the multiplicative noise problem using the Stratonovich or a
Stratonovich-like convention. The Itô-convention clearly does not give the
Boltzmann distribution for equilibrium and it gives the MBD. The dilemma of
Itô-vs-Stratonovich conventions is apparently resolved about 40 years ago by
the understanding of the fact that different conventions represent different
physical processes [20, 21]. Therefore, the sole purpose in the context of
thermal equilibrium should be to identify which convention corresponds to a
process which is consistent with the demands of a thermal equilibrium. One may
also note, in this context, that the Stratonovich or Stratonovich-like
conventions involve correlated (anticipating) noise. There have also been
attempts to look into this issue in some alternative way [22] where it is
argued that specification of coordinate dependence of diffusivity needs to be
accompanied by additional information to specify the state. The Itô,
Stratonovich and Hänggi-Klimontovich conventions have been related to an
infinite density and its shape by Leibovich and Barkai in [23] to shed light
on various interpretations.

In the present paper, we focus on the results of a classic experiment reported
in 1994 by Faucheux and Libchaber [24] to compare the BD and the MBD with the
results of this experiment. We particularly choose this experimental result
for comparison of the BD with the MBD because, one can do this comparison here
without involving any free (adjustable) parameter for curve fitting. When two
different functions are compared with experimental outcome, ideally, one must
not use any free parameter for curve fitting and this experimental result of
[24] gives us that opportunity. The comparison shows good matching of the
experimental results of average coordinate-dependent diffusivity with those
computed using the MBD.

To our knowledge, it is the first experiment that had attempted to find the
trajectory of a BP by direct visualization of the particle [24]. Based on this
direct tracking, the authors have made an effort to reconcile the position-
dependent-diffusivity to the Stokes-Einstein relation i.e., they have
considered thermal equilibration of the BP. Although, the BD had been used,
however, there are observations reported in the same paper that the BD could
not fit experimental data in the way it should have fitted it. Moreover, each
and every parameter of the system is experimentally determined and reported in
the paper for each and every particle tracked [24].

The plan of the present paper goes in the following way. In the next section,
we would first describe the experiment that has been done by Faucheux and
Libchaber and we give a description of the theoretical methods used by these
authors. In the next section we will describe the details of comparison of the
experimental results with the BD and the MBD where the theory practically
remains the same as the one used by Faucheux and Libchaber except for the
method of their theoretical curve fitting. Following that, we conclude the
paper with a discussion.

##  II The experiment and theoretical method of Faucheux and Libchaber

###  II.1 The experiment

In the experiment of Faucheux and Libchaber [24], they observed Brownian
motion of silica and latex beads (considered spheres) of diameter 1 to 3
$\mu$m in ultra-pure water confined within two clean horizontal glass plates
of spacing between 6 to 1000 $\mu$m. The real-time position of the diffusing
spheres on a horizontal plane ($x$ and $y$ coordinate) were taken by a CCD
camera coupled to the microscope. A time series of about 3 minute were
constructed for each of these beads where the typical diffusion coefficient is
about 1 $\mu$m2/s. Typically, 12 such time series for each bead have been
considered for an ensemble average which corresponds to about half an hour of
recording for each bead. The typical time over which the water gets
contaminated in this experiment was assessed to be about a day.

Based on these time series of the horizontal positions of individual beads the
diffusivity $D_{||}$ of the particles were determined over a horizontal plane
which actually is a projection of the 3-dimensional trajectory of the bead on
a plane perpendicular to the vertical. The diffusivity on the horizontal plane
$D_{||}$ and the diffusivity of the bead in the vertical ($z$-direction)
$D_{\perp}$ are functions of the distance of the bead from the confining
horizontal glass plates. Therefore, the horizontal diffusivity of a bead
measured in this projection is an average over the vertical direction
($z$-coordinate) where the probability density for position distribution could
be the BD or the MBD which is to be found out. Care was taken in this
experiment to avoid the beads coming too close to each other and interact
hydro-dynamically such that the diffusion remains only a function of the
vertical $z$-coordinate. For other details of the experimental determination
of the $D_{||}$ refer Faucheux and Libchaber [24].

Figure 1: Experimentally determined ratio $D_{||}/D_{0}$ as a function of
$\gamma$

###  II.2 Coordinate dependence of horizontal-diffusivity on the distance from
wall

Among these beads, the heavier ones would have a smaller average distance from
the lower glass plate where the lighter ones would have average height
relatively larger. The diffusive excursion of these beads will be around these
average height in such a way that the average height remains stationary under
steady-state (thermal equilibrium) conditions. Keeping this in mind, Faucheux
and Libchaber used a parameter defined as

$\gamma=\frac{h-r}{r}$ (2)

where $r$ is the radius of the bead and the average vertical height $h$ is
defined as

$\displaystyle h$ $\displaystyle=$ $\displaystyle\int_{r}^{t-r}{zP_{B}(z)dz}$
(3) $\displaystyle=$
$\displaystyle\frac{e^{-r/L}(rL+L^{2})-e^{(r-t)/L}([t-r]L+L^{2})}{L(e^{-r/L}-e^{(r-t)/L})},$

where

$P_{B}(z)=\left(\frac{1}{L}\right)\frac{e^{-z/L}}{e^{-r/L}-e^{(r-t)/L}},$ (4)
$L=k_{B}T/g\bigtriangleup\\!m$ (5)

and

$\bigtriangleup m=\frac{4}{3}\pi r^{3}(\rho-\rho_{0}),$ (6)

to plot the experimentally obtained $D_{||}$ against this parameter. This is a
nice choice of a parameter for the horizontal diffusivity $D_{||}$ to be
plotted against because there is vertical excursion of the bead associated to
its horizontal motion and its vertical as well as horizontal diffusivity are
functions of the vertical distance of the bead from the confining plates.
However, in the experimentally evaluated $D_{||}$ this dependence of
$z$-coordinate appears on an average in this $z$-direction because $D_{||}$ is
evaluated from a projection on the $x-y$ plane. In the above expressions
$\rho$ is the density of the material of the bead, $\rho_{0}$ is the density
of water considered in the present paper to be $1000$ kg/m3 at room
temperature $T=300$ K, $g$ is the acceleration due to gravity taken to be
$9.8$ m/s2, $t$ is the spacing between the confining plates and $z$ is the
vertical coordinate of the centre of mass of the bead taken to be a sphere.

The experimentally measured $D_{||}$ normalized by the bulk diffusivity
$D_{0}$ is then plotted against this parameter $\gamma$ (following [24]) in
Fig.1. In this figure, the data-points (as given in the Table I of [24]) are
plotted by the symbol of square boxes and the error-bars are extracted from
the Fig.4 of [24] using WebPlotDigitizer (https://apps.automeris.io/wpd/). In
what follows, the same box symbol and error-bars would be used in all the
figures to correspond to the experimental data-points. Relevant data from the
Table I of the paper [24] are shown in the Table I of the present paper.

Sample No. |  Bead diameter ($\mu$m) |  Bead density ($g/cm^{3}$) |  L ($\mu$m) |  Sample thickness t ($\mu$m) |  $D_{||}/D_{0}$ (expt) |  $\gamma$
---|---|---|---|---|---|---
1 | 2.5 | 2.1 | 0.05 | 6 | 0.32 | 0.04
2 | 2.5 | 2.1 | 0.05 | 25 | 0.36 | 0.04
3 | 2.5 | 2.1 | 0.05 | 100 | 0.36 | 0.04
4 | 2.5 | 2.1 | 0.05 | 1000 | 0.36 | 0.04
5 | 3.5 | 1.05 | 0.37 | 12 | 0.4 | 0.2
6 | 3.5 | 1.05 | 0.37 | 50 | 0.41 | 0.2
7 | 3 | 1.05 | 0.6 | 12 | 0.35 | 0.4
8 | 2.5 | 1.05 | 1.03 | 12 | 0.38 | 0.8
9 | 2 | 1.05 | 2.02 | 6 | 0.45 | 1.4
10 | 1 | 2.1 | 0.73 | 6 | 0.41 | 1.4
11 | 1 | 2.1 | 0.73 | 25 | 0.42 | 1.4
12 | 1 | 2.1 | 0.73 | 100 | 0.56 | 1.4
13 | 2 | 1.05 | 2.02 | 12 | 0.61 | 1.9
14 | 1.5 | 1.05 | 4.8 | 6 | 0.52 | 2.5
15 | 1.5 | 1.05 | 4.8 | 12 | 0.78 | 4.6
16 | 1 | 1.05 | 16.1 | 6 | 0.71 | 4.6
17 | 1.5 | 1.05 | 4.78 | 25 | 0.81 | 6.1
18 | 1 | 1.05 | 16.1 | 12 | 0.87 | 9.4
19 | 1 | 1.05 | 16.1 | 25 | 0.82 | 17
20 | 1 | 1.05 | 16.1 | 50 | 0.89 | 23

Table 1: This table contains the relevant data as is given in the Table 1 of
[24]. This data is used to plot $D_{||}/D_{0}$ vs $\gamma$ in Fig.1.

###  II.3 Effect of the wall on the damping and diffusivity

The effect of the wall on the damping coefficients (as considered by Faucheux
and Libchaber [25, 26, 27, 28]) are well known with support from molecular-
dynamics [29].

$\displaystyle\eta_{x}$ $\displaystyle=\eta_{y}\simeq$ $\displaystyle\eta_{0}$
$\displaystyle[1-\frac{9}{16}(r/z)+\frac{1}{8}(r/z)^{3}-\frac{45}{256}(r/z)^{4}-\frac{1}{16}(r/z)^{5}]^{-1},$

and

$\displaystyle\eta_{z}$ $\displaystyle=$
$\displaystyle\frac{4}{3}\eta_{0}\sinh{\alpha}\sum_{n=1}^{\infty}\frac{n(n+1)}{(2n-1)(2n+3)}\times$
$\displaystyle\left[\frac{2\sinh{(2n+1)}\alpha+(2n+1)\sinh{2\alpha}}{4\sinh^{2}{(n+1/2)\alpha}-(2n+1)^{2}\sinh^{2}{\alpha}}-1\right],$

where $\alpha=\cosh^{-1}(z/r)$ and $\eta_{0}$ is the damping coefficient in
the bulk of the fluid which we will consider in this paper to have a value
$\eta_{0}=8.50\times 10^{-4}$ Pa.s. The relation (7) is exact and relation (6)
is accurate up to the order $(r/z)^{5}$. In this paper, we present results of
computations up to $n=5$ while using (7) to keep parity with the order of
accuracy with (6), however, we have checked that going to a few higher orders
in $n$ does not change the result appreciably.

(a)

(b)

Figure 2: Theoretically calculated $\langle D_{||}\rangle^{B/MB}/D_{0}$
compared with experimental data points using $\gamma_{B/MB}$ as parameters.
Theoretical data points are fitted using cubic splines with smoothing
parameter as 0.01

Theoretical diffusivity on the horizontal plane is thus considered to be
$D_{h}=k_{B}T/6\pi\eta_{x}r=k_{B}T/6\pi\eta_{y}r$ and that in the vertical
direction is $D_{\perp}=k_{B}T/6\pi\eta_{z}r$, where $D_{h}=D_{x}=D_{y}$ and
$D_{\perp}=D_{z}$. These expressions of the theoretical diffusivities follow
from the standard fluctuation-dissipation relation or the Stokes-Einstein
relation under thermal equilibrium conditions. In this particular experiment,
the vertical separation ($t$) between the confining horizontal glass plates
vary from 6 $\mu$m to 1000 $\mu$m. It requires to take into account the
presence of both the plates and the vertical distance of the bead from them
when the separation $t$ is small. Taking into account the effect of both the
plates at the bottom and at the top (upper plate), Faucheux and Libchaber [24]
introduced an effective $\eta_{x,y,z}$ as

$\eta_{x,y,z}=\eta^{b}_{x,y,z}+\eta^{u}_{x,y,z}-\eta_{0},$ (9)

where the superscripts $b$ and $u$ indicate damping coefficients due to bottom
and upper (top) plates respectively. Note that, $\eta^{b}_{x,y,z}$ is a
function of $z/r$ where $\eta^{u}_{x,y,z}$ is similarly a function of
$(t-z)/r$. This expression of damping coefficient $\eta_{x,y,z}$ has nice
limits. When the upper plate goes to infinity $(t\to\infty)$, we get
$\eta^{u}_{x,y,z}=\eta_{0}$ and only the lower plate has any effect on the
bead. Thus, even when the bottom and the top plates are far away from the BP
we get $\eta_{x,y,z}=\eta_{0}$.

The relation (8) takes into account the independent presence of the bottom and
the top plate on the damping coefficient of the BP without considering any
interaction. Imagine that, only the bottom plate is present and the damping
coefficient is modified due only to its presence such that

$\eta_{x,y,z}=\eta_{0}+\bigtriangleup\\!\eta^{b}(z),$ (10)

where $\bigtriangleup\\!\eta^{b}(z)$ is the modification of the damping due to
the presence of the bottom plate on top of what it is there in the bulk.
Similarly, when one takes into account the presence of the top plate whose
hydrodynamic effect is not going to interfere with that produced by the bottom
plate, there would also be an increase in the damping by the amount
$\bigtriangleup\\!\eta^{u}(z)$ on top of what the bottom plate has produced.
Therefore the new damping coefficient considering the presence of both the
plates is

$\displaystyle\eta_{x,y,z}$ $\displaystyle=$
$\displaystyle\eta_{0}+\bigtriangleup\\!\eta^{b}(z)+\bigtriangleup\\!\eta^{u}(z)$
$\displaystyle=$
$\displaystyle[\eta_{0}+\bigtriangleup\\!\eta^{b}(z)]+[\eta_{0}+\bigtriangleup\\!\eta^{u}(z)]-\eta_{0}$
$\displaystyle=$ $\displaystyle\eta^{b}_{x,y,z}+\eta^{u}_{x,y,z}-\eta_{0}$

We have found this expression (8) to work quite well and in what follows we
would go by the relation (8) to take into account the hydrodynamic effects due
to both the walls.

##  III Comparison of BD and MBD with the experimental result

In this section we compare the analytically obtained height-averaged
horizontal diffusivity $\langle D_{||}\rangle^{B/MB}$ with the experimentally
obtained results ($D_{||}$) of Faucheux and Libchaber [24]. We use two
different distributions the BD and the MBD to get the height-average
($z$-averaged) horizontal diffusivity in the following way

$\displaystyle\langle D_{||}\rangle^{B}$ $\displaystyle=$
$\displaystyle\int_{r+\epsilon}^{t-(r+\epsilon)}{D_{h}(z)P_{B}(z)dz}$ (11)
$\displaystyle\langle D_{||}\rangle^{MB}$ $\displaystyle=$
$\displaystyle\int_{r+\epsilon}^{t-(r+\epsilon)}{D_{h}(z)P_{MB}(z)dz},$ (12)

The distributions are

$\displaystyle P_{B}$ $\displaystyle=$ $\displaystyle N_{B}\exp{(-z/L)}$ (13)
$\displaystyle P_{MB}$ $\displaystyle=$
$\displaystyle\frac{N_{MB}}{D_{\perp}(z)}\exp{(-z/L)},$ (14)

where $N_{B}$ and $N_{MB}$ are the normalization constants for the BD and the
MBD respectively, evaluated exactly over the same lower/upper limits for a
particular bead as shown in (9) and (10). We implement a cut-off distance
$\epsilon=10^{-3}\mu$m to avoid the divergence appearing in the MBD due to the
vanishing of $D_{\perp}(z)$ at $z=r$ and $z=t-r$. We implement the same cut-
off in (9) and (10) in order to make the averages with respect to the BD and
the MBD identically computed. The integration here is a left-Riemann-sum with
bin-width $10^{-3}\mu$m which is being presented. We have checked that the
left and right Riemann-sum are giving excellent match with this bin-size and
this is to be checked because of the divergence near the walls.

Sample No. |  $h_{B}$ ($\mu$m) |  $h_{MB}$ ($\mu$m) |  $\gamma_{B}$ |  $\gamma_{MB}$ |  $\Gamma_{MB}$ |  $tL/100$ $\mu$ $m^{2}$ |  $\langle D_{||}\rangle^{B}/D_{0}$ |  $\langle D_{||}\rangle^{MB}/D_{0}$ |  $D_{||}/D_{0}$ (expt)
---|---|---|---|---|---|---|---|---|---
1 | 1.30 | 1.26 | 0.04 | 0.01 | 0.01 | 0.003 | 0.34 | 0.32 | 0.32
2 | 1.30 | 1.26 | 0.04 | 0.01 | 0.01 | 0.0125 | 0.36 | 0.33 | 0.36
3 | 1.30 | 1.26 | 0.04 | 0.01 | 0.01 | 0.05 | 0.37 | 0.34 | 0.36
4 | 1.30 | 1.26 | 0.04 | 0.01 | 0.01 | 0.5 | 0.37 | 0.34 | 0.36
5 | 2.12 | 1.83 | 0.21 | 0.04 | 0.04 | 0.0444 | 0.45 | 0.35 | 0.4
6 | 2.12 | 1.83 | 0.21 | 0.04 | 0.04 | 0.185 | 0.47 | 0.36 | 0.41
7 | 2.10 | 1.64 | 0.40 | 0.09 | 0.09 | 0.072 | 0.51 | 0.37 | 0.35
8 | 2.28 | 1.52 | 0.82 | 0.21 | 0.21 | 0.1236 | 0.59 | 0.41 | 0.38
9 | 2.38 | 1.80 | 1.38 | 0.80 | 0.69 | 0.1212 | 0.60 | 0.45 | 0.45
10 | 1.23 | 0.75 | 1.45 | 0.49 | 0.47 | 0.0438 | 0.66 | 0.47 | 0.41
11 | 1.23 | 0.73 | 1.46 | 0.46 | 0.45 | 0.1825 | 0.69 | 0.48 | 0.42
12 | 1.23 | 0.73 | 1.46 | 0.46 | 0.45 | 0.73 | 0.69 | 0.48 | 0.56
13 | 2.95 | 1.73 | 1.95 | 0.73 | 0.68 | 0.2424 | 0.69 | 0.49 | 0.61
14 | 2.65 | 2.23 | 2.54 | 1.97 | 1.48 | 0.288 | 0.67 | 0.52 | 0.52
15 | 4.22 | 2.99 | 4.63 | 2.99 | 2.43 | 0.576 | 0.78 | 0.60 | 0.78
16 | 2.87 | 2.64 | 4.74 | 4.28 | 2.75 | 0.966 | 0.75 | 0.60 | 0.71
17 | 5.36 | 3.18 | 6.14 | 3.24 | 2.93 | 1.195 | 0.83 | 0.63 | 0.81
18 | 5.38 | 4.82 | 9.76 | 8.63 | 5.53 | 1.932 | 0.84 | 0.71 | 0.87
19 | 9.62 | 8.28 | 18.25 | 15.56 | 10.72 | 4.025 | 0.90 | 0.80 | 0.82
20 | 14.15 | 11.64 | 27.30 | 22.29 | 17.32 | 8.05 | 0.93 | 0.84 | 0.89

Table 2: Contains values of different parameters ($\gamma_{B}$, $\gamma_{MB}$,
$\Gamma_{MB}$ and $tL/100$) and diffusivity ratios, both theoretical ($\langle
D_{||}\rangle^{B/MB}/D_{0}$) and experimental $D_{||}/D_{0}$. Also, $h_{B}$
and $h_{MB}$ are provided as they are used to calculate the parameters.

It is important to note that, for a matching of the experimental results with
the theoretical ones (9) and (10), we need to find out the theoretical values
of the average diffusion coefficients for each individual particle. The
corresponding parameter $\gamma$ has also to be found out for each individual
particle because the boundary conditions, matter density and size are quite
different for different particle which actually correspond to different
experimental conditions. In Fig.2(a) and (b), we have compared the theoretical
values (according to (9) and (10)) of the average diffusivity of individual
particles against the experimental data using the parameter $\gamma_{B}$ and
$\gamma_{MB}$. These are the same $\gamma$ parameter of Faucheux and
Libchaber, but, computed by us within the specified limits using BD and MBD
respectively. Theoretical data points for Fig.2 and Fig.3 are given in Table
II. The continuous line in each figure corresponds to the cubic-spline to the
theoretical data points with a smoothing factor $0.01$. As we can see that,
the parameters $\gamma_{B/MB}$ are distribution dependent, and that is why we
have plotted two different graphs and have not overlaid the theoretical
diffusivities on the same graph.

Faucheux and Libchaber have pointed out in ref.[24] that the BD does not work
to match the experimental data when used as in (9). Fig.2(a) shows that this
assessment of Faucheux and Libchaber is correct and the Boltzmann distribution
fails to match the trend of the experimental result particularly in the region
of lower values of $D_{||}$ where the modification of the diffusivity is
appreciable due to the proximity of the wall. The theoretical and experimental
values are quite close while the BD is used at the two ends i.e., for very
heavy particles which practically stay at the bottom and the ones which are
very light and are mostly in the bulk. Note that, very heavy beads which stay
near the bottom plate they do not in general undergo much vertical excursion
and as a result the effect of averaging in the vertical direction is quite
small for these particles. Moreover, very light ones are staying mostly in the
bulk of the fluid where the diffusivity is less perturbed and thereby the
effect of average in the vertical direction would have little effect on those
as well. It is the trend in the intermediate region which is the most
important part. However, in this middle region, the trend of the variation of
the experimental result is missed by $\langle D_{||}\rangle^{B}$.

On the contrary, Fig.2(b) shows a much better match between the theoretical
and experimental result in the intermediate region when the MBD is used. It
also appears that the MBD is not matching the larger diffusivity region as
much as the BD does. Note that, the error bars in the larger diffusivity
region are larger due to limitations in the averaging process over the
experimental trajectories mentioned in [24]. Let us look at the data of the
larger diffusivity (near bulk) more closely. For example, the sample number 18
in Table I, which is the third experimental data-point in the plots from the
right hand side end is clearly an over-count. This particle has a bead
diameter 1 $\mu$m which is the same as the bead of the sample number 19 and
both have the same density. Where the sample number 18 has a sample thickness
$t=12$ $\mu$m that of the sample number 19 is 25. Given this, the
$D_{||}/D_{0}$ of sample number 18 should be smaller than that of the sample
number 19, however, that is not the experimental result as shown in Table I.
Similarly, if one compares the experimental value of $D_{||}/D_{0}$ of sample
number 17 with that of the sample number 19, one would find that the
experimental value is larger (quite close to that of sample number 19) than
that it should be for the sample number 17 given the bigger diameter of the
bead. The experimental data in Table I, in this way, reveals that there
apparently exists some bias towards over-counting the diffusivity of the
lighter particles.

(a)

(b)

Figure 3: Theoretical $\langle D_{||}\rangle^{MB}/D_{0}$ compared with
experimental data points and $\langle D_{||}\rangle^{B}/D_{0}$ using
parameters $\Gamma_{MB}$ and $tL/100$, respectively, in (a) the cubic spline
to the theoretical result has a smoothing factor 0.01 and in (b) A cubic
spline with smoothing factor 0.1 is drawn for $\langle
D_{||}\rangle^{MB}/D_{0}$.

In the Fig.3(a) and (b), we compare the experimental results of Faucheux and
Libchaber [24] with the $\langle D_{||}\rangle^{MB}$ where we use two
different parameters $\Gamma_{MB}$ and $tL/100$, respectively. The latter
parameter is a distribution independent parameter and we can overlay graphs
corresponding to the BD and the MBD here. We have used the modification of the
Faucheux and Libchaber’s parameter as $\Gamma_{MB}$ (in Fig.3(a)) to take into
account the distance of the particle from both the lower and the upper
confining planes in the parameter as

$\Gamma_{MB}=\gamma_{MB}\gamma^{\prime}_{MB},$ (15)

where

$\gamma^{\prime}_{MB}=1-\gamma_{MB}\frac{r}{t},$ (16)

where $h$ is evaluated by averaging with MBD. $\gamma^{\prime}_{MB}$ is a
fraction which has the limit of unity as the upper plate goes far away
$(t\to\infty)$ and the effect of the lower plate only is there on the
parameter at that limit. Otherwise, the presence of the factor
$\gamma^{\prime}_{MB}$ lowers the value of the parameter $\Gamma_{MB}$ in
comparison with $\gamma_{MB}$ depending upon the proximity of the upper plate.
As has already been pointed out by Faucheux and Libchaber, effect of the upper
plate in the parameter $\gamma$ is not felt that much is clear from the
comparison of the Fig.2(b) and Fig.3(a).

The parameter in Fig.3(b) is $tL/100$ (normalization is done to have similar
scales in all the figures) which is distribution independent. In this figure
we have plotted $\langle D_{||}\rangle^{B}$ by green dots and the $\langle
D_{||}\rangle^{MB}$ by red stars. We have drawn the cubic spline to only the
$\langle D_{||}\rangle^{MB}$ and not for $\langle D_{||}\rangle^{B}$ because
the data points are too scattered. This figure clearly shows us that the
matching of the theory with experimental data remains quite good with the MBD
for the beads staying closer to the lower plate (and hence feeling more
hydrodynamic effects) in comparison to those which are lighter. However, as we
have already discussed above, from the table of experimental data points for
$D_{||}$ as is given by Faucheux and Libchaber, there probably exists a trend
of a little over-counting for the $D_{||}$ of the lighter particles and this
could be resolved by a similar experiment if done at present times.

Corresponding to the parameter $tL/100$ the sample number 2, 3 and 4 sit at
the parameter values $tl/100=0.0125$, $0.05$ and $0.5$. However, all these
samples being heaviest particles have the same experimentally measured
$D_{||}/D_{0}$ indicating that the sample thickness $t$ beyond $t=25$ $\mu$m
does not have any effect on the diffusivity of these particles because those
are anyway very close to the bottom wall. However, the parameter $tL/100$ is
making these sample points sit at different points on the horizontal axis
because the parameter counts the sample thickness even where it is not
influencing the result. Apart from this, the parameter $tL/100$ appears to be
a good parameter for almost all other sample numbers to directly compare the
effects of the BD and the MBD with the experimental results.

##  IV Discussion

The thermal equilibrium distribution of a BP with coordinate dependent
diffusion is the MBD has been theoretically argued for in the references [11,
12]. In the reference [13] it has been shown that this coordinate-dependent
diffusion is indeed a continuum limit of a random walk in contact with a heat-
bath. The basis of the theoretical argument in favour of an Itô-process as
representing thermal equilibrium has been identifying the Itô-process to be
consistent with the local applicability of the Stokes-Einstein relation
(fluctuation-dissipation relation). Moreover, an Itô-process does not
introduce any correlation in the thermal noise which is unavoidable when the
Stratonovich or a Stratonovich-like conventions are used.

A comparison of the Itô-process with existing experimental result can really
set the stage for further serious exploration in this direction. There are
number of ways in which the BD is theoretically imposed for equilibrium using
the Stratonovich or a Stratonovich-like convention in such a multiplicative
noise problems. These theories are inconsistent (despite imposing the BD) from
the stand point of fundamental requirements of a stochastic process to be a
thermal equilibrium. We have given here an alternative explanation to the
experimental results of [24] which had an explanation already given by
Faucheux and Libchaber using the BD. Moreover, it is apparent from this
comparison that the MBD (corresponding to an Itô-process) better captures the
experimental results than the BD.

It is crucial for fundamental reasons that the equilibrium distribution of a
BP with coordinate-dependent (or state-dependent for a many-body system)
diffusion be experimentally settled. Coordinate- or state-dependent diffusion
has the potential to generalize our knowledge of thermal equilibrium to a
great extent because this process takes into account an in-homogeneity of
space keeping intact the thermal homogeneity of the heat-bath. In the
conventional theory of uniform diffusion, this in-homogeneity of space is
absent. In the presence of coordinate-dependent diffusion the homogeneity of
space (even in the absence of a conservative force) gets broken and decouples
from the thermal homogeneity of the heat-bath so that the particle can see a
constant temperature in an otherwise in-homogeneous space. This decoupling can
result in counter-intuitive reality to our present understanding of thermal
equilibrium.

Diffusion is a major source of transport in numerous biological systems and
the space remains quite crowded in such cases where the coordinate dependence
of the diffusion could be the rule and not an exception. Moreover, coordinate
dependence of diffusion could have other counter-intuitive effects of
producing directed motions in structured objects as is shown in [30, 31].
These two effects of coordinate dependence of diffusion i.e. modifying the
Boltzmann distribution and spontaneously filtering of Brownian motion by
structured objects in the presence of state dependent diffusion may prove to
be essential ingredients for our understanding of a plethora of other
phenomena that happen at meso-scales being driven by equilibrium fluctuations.

The experimental results of [24] is very reliable because it has used ultra-
pure water as the fluid. The use of ultra-pure water and other experimental
cares taken to keep the fluid pure (or doing the experiment within the time
the fluid remains contaminant-free) has kept the potential term for the
Boltzmann factor to be in its simplest form. Practically no amount of free
ions were presumably present in the fluid helping avoid the the presence of
any Debye-Huckel term. This is the reason that there is no need to deal with
any free parameter in this context. The experiment was done about 30 years
ago, however, still it remains ever more important because of the way it was
done. At present, the particle tracking ability has improved to an extremely
high level of accuracy. It is a matter of time, to our belief, that the MBD
could be shown to be the equilibrium distribution for coordinate-dependent
diffusion as opposed to the BD once one does a similar experiment where there
would be no need to use adjustable parameters to fit the experimental results.

## References

  * Roussel and Roussel [2004] C. J. Roussel and M. R. Roussel, Reaction–diffusion models of development with state-dependent chemical diffusion coefficients, Progress in biophysics and molecular biology 86, 113 (2004).
  * Barik and Ray [2005] D. Barik and D. S. Ray, Quantum state-dependent diffusion and multiplicative noise: a microscopic approach, Journal of Statistical Physics 120, 339 (2005).
  * Sargsyan  _et al._ [2007] V. Sargsyan, Y. V. Palchikov, Z. Kanokov, G. Adamian, and N. Antonenko, Coordinate-dependent diffusion coefficients: Decay rate in open quantum systems, Physical Review A 75, 062115 (2007).
  * Chahine  _et al._ [2007] J. Chahine, R. J. Oliveira, V. B. Leite, and J. Wang, Configuration-dependent diffusion can shift the kinetic transition state and barrier height of protein folding, Proceedings of the National Academy of Sciences 104, 14646 (2007).
  * Best and Hummer [2010] R. B. Best and G. Hummer, Coordinate-dependent diffusion in protein folding, Proceedings of the National Academy of Sciences 107, 1088 (2010).
  * Lai  _et al._ [2014] Z. Lai, K. Zhang, and J. Wang, Exploring multi-dimensional coordinate-dependent diffusion dynamics on the energy landscape of protein conformation change, Physical Chemistry Chemical Physics 16, 6486 (2014).
  * Berezhkovskii and Makarov [2017] A. M. Berezhkovskii and D. E. Makarov, Communication: Coordinate-dependent diffusivity from single molecule trajectories, The Journal of chemical physics 147, 201102 (2017).
  * Foster  _et al._ [2018] D. A. Foster, R. Petrosyan, A. G. Pyo, A. Hoffmann, F. Wang, and M. T. Woodside, Probing position-dependent diffusion in folding reactions using single-molecule force spectroscopy, Biophysical journal 114, 1657 (2018).
  * Ghysels  _et al._ [2017] A. Ghysels, R. M. Venable, R. W. Pastor, and G. Hummer, Position-dependent diffusion tensors in anisotropic media from simulation: oxygen transport in and through membranes, Journal of chemical theory and computation 13, 2962 (2017).
  * Yamilov  _et al._ [2014] A. G. Yamilov, R. Sarma, B. Redding, B. Payne, H. Noh, and H. Cao, Position-dependent diffusion of light in disordered waveguides, Physical review letters 112, 023904 (2014).
  * Bhattacharyay [2019] A. Bhattacharyay, Equilibrium of a brownian particle with coordinate dependent diffusivity and damping: Generalized boltzmann distribution, Physica A: Statistical Mechanics and its Applications 515, 665 (2019).
  * Bhattacharyay [2020] A. Bhattacharyay, Generalization of stokes–einstein relation to coordinate dependent damping and diffusivity: an apparent conflict, Journal of Physics A: Mathematical and Theoretical 53, 075002 (2020).
  * Maniar and Bhattacharyay [2021] R. Maniar and A. Bhattacharyay, Random walk model for coordinate-dependent diffusion in a force field, Physica A: Statistical Mechanics and its Applications 584, 126348 (2021).
  * Sokolov [2010] I. M. Sokolov, Itô, stratonovich, hänggi and all the rest: The thermodynamics of interpretation, Chemical Physics 375, 359 (2010).
  * Lau and Lubensky [2007] A. W. Lau and T. C. Lubensky, State-dependent diffusion: Thermodynamic consistency and its path integral formulation, Physical Review E 76, 011123 (2007).
  * Sancho  _et al._ [1982] J. Sancho, M. S. Miguel, and D. Dürr, Adiabatic elimination for systems of brownian particles with nonconstant damping coefficients, Journal of Statistical Physics 28, 291 (1982).
  * Sancho [2011] J. M. Sancho, Brownian colloidal particles: Ito, stratonovich, or a different stochastic interpretation, Physical Review E 84, 062102 (2011).
  * Farago and Grønbech-Jensen [2014a] O. Farago and N. Grønbech-Jensen, Langevin dynamics in inhomogeneous media: Re-examining the itô-stratonovich dilemma, Physical Review E 89, 013301 (2014a).
  * Farago and Grønbech-Jensen [2014b] O. Farago and N. Grønbech-Jensen, Fluctuation–dissipation relation for systems with spatially varying friction, Journal of Statistical Physics 156, 1093 (2014b).
  * Van Kampen [1981] N. G. Van Kampen, Itô versus stratonovich, Journal of Statistical Physics 24, 175 (1981).
  * Mannella and McClintock [2022] R. Mannella and P. V. McClintock, Itô versus stratonovich: 30 years later, The Random and Fluctuating World: Celebrating Two Decades of Fluctuation and Noise Letters , 9 (2022).
  * Tupper and Yang [2012] P. Tupper and X. Yang, A paradox of state-dependent diffusion and how to resolve it, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 468, 3864 (2012).
  * Leibovich and Barkai [2019] N. Leibovich and E. Barkai, Infinite ergodic theory for heterogeneous diffusion processes, Physical Review E 99, 042138 (2019).
  * Faucheux and Libchaber [1994] L. P. Faucheux and A. J. Libchaber, Confined brownian motion, Physical Review E 49, 5158 (1994).
  * Faxen [1924] H. Faxen, Ark. Mat. Astron. Fys. 18, 1 (1924).
  * Brenner [1961] H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, Chemical Engineering Science 16, 242 (1961).
  * Cox and Brenner [1967] R. G. Cox and H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface. ii - small gap widths, including inertial effects., Chemical Engineering Science 22, 1753 (1967).
  * Goldman  _et al._ [1967] A. J. Goldman, R. G. Cox, and H. Brenner, Slow viscous motion of a sphere parallel to a plane wall—i motion through a quiescent fluid, Chemical engineering science 22, 637 (1967).
  * Chio and Tse [2020] C. C. Chio and Y.-L. S. Tse, Hindered diffusion near fluid–solid interfaces: Comparison of molecular dynamics to continuum hydrodynamics, Langmuir 36, 9412 (2020).
  * Bhattacharyay [2012] A. Bhattacharyay, Directed transport in equilibrium: A model study, Physica A: Statistical Mechanics and its Applications 391, 1111 (2012).
  * Sharma and Bhattacharyay [2020] M. Sharma and A. Bhattacharyay, Conversion of heat to work: An efficient inchworm, Physica Scripta 95, 105004 (2020).