Large Deviations in One-Dimensional Random Sequential Adsorption

Here we outline the computation of the average number $A_{L}=\langle N\rangle$ of occupied sites on an open interval of length $L$. We employ an approach essentially developed by Flory F39 , see also a textbook book ; our analysis (Sec. III) of the cumulant generating function will proceed along similar lines.

For small $L$, one can compute $A_{L}=\langle N\rangle$ by hand. For instance, $N=1$ when $L=1$ and 2; for $L=4$, one gets $N=2$. For $L=3$ and $L\geq 5$, the total number $N$ of absorbed particles varies from realization to realization, e.g. $N=1$ (probability $1/3$) or $N=2$ (probability $2/3$) when $L=3$. Thus

$$A_{1}=1,\quad A_{2}=1,\quad A_{3}=\frac{5}{3},\quad A_{4}=2$$

(2)

To compute $A_{L}$ for arbitrary $L$, suppose $k$ is the landing site of the first particle. The intervals of lengths $k-2$ and $L-k-1$ on the left and right of the first deposited particle are subsequently filled independently. This key feature makes the one-dimensional situation tractable. For our model, we arrive at the recurrence

$$A_{L}=\frac{1}{L}\sum_{k=1}^{L}\big{(}A_{k-2}+1+A_{L-k-1}\big{)}$$

(3)

that is applicable for all $L\geq 1$ if we set $A_{-1}=A_{0}=0$. Using the generating function

$$A(x)=\sum_{L\geq 1}A_{L}\,x^{L}$$

(4)

we convert the recurrence (3) into a differential equation

$$\frac{dA(x)}{dx}=\frac{2x}{1-x}\,A(x)+\frac{1}{(1-x)^{2}}$$

(5)

Solving this linear inhomogeneous differential equation subject to the initial condition $A(0)=0$, we obtain

$$A(x)=\frac{1-e^{-2x}}{2(1-x)^{2}}$$

(6)

from which

$$A_{L}=\sum_{k=0}^{L-1}\frac{(-2)^{k}}{(k+1)!}\,(L-k)$$

(7)

Using (7) one can extract the large $L$ behavior:

$$A_{L}=\frac{1-e^{-2}}{2}\,(L+3)-1+\frac{(-2)^{L+1}}{(L+2)!}+\ldots$$

(8)

The average jamming density is $\rho_{L}=A_{L}/L$, and in the $L\to\infty$ limit the jamming density becomes

$$\rho_{\text{jam}}=\rho_{\infty}=\frac{1-e^{-2}}{2}=0.432332358\ldots$$

(9)

When $L\gg 1$, the probability that a site in the bulk is occupied is very close to (9). Near the boundary the densities are different.

Let us first compute the probability $p_{L}$ that the leftmost site is occupied. For small $L$, one can compute $p_{L}$ by hand to yield

$$p_{1}=1,\quad p_{2}=\frac{1}{2},\quad p_{3}=\frac{2}{3},\quad p_{4}=\frac{5}{8}$$

(10)

etc. Generally, the probability $p_{L}$ satisfies the recurrence

$$p_{L}=\frac{1}{L}+\frac{1}{L}\sum_{k=1}^{L-2}p_{k}$$

(11)

Indeed, with probability $L^{-1}$, the first landing site is the leftmost site; this explains the first term on the right-hand side of (11). The leftmost site may also get occupied if the first landing site is $k+1$ with $k=1,\ldots,L-2$. Each such event happens with probability $L^{-1}$, and the leftmost site in the segment of length $k$ could be occupied with probability $p_{k}$. This explains the sum on the right-hand side of (11).

Using the generating function

$$P(x)=\sum_{L\geq 1}p_{L}\,x^{L}$$

(12)

we convert the recurrence (11) into a differential equation

$$\frac{dP(x)}{dx}=\frac{x}{1-x}\,P(x)+\frac{1}{1-x}$$

(13)

which is solved to yield

$$P(x)=\frac{1-e^{-x}}{1-x}$$

(14)

Expanding $P(x)$ we obtain

$$p_{L}=\sum_{k=1}^{L}\frac{(-1)^{k-1}}{k!}$$

(15)

In the $L\to\infty$ limit, i.e. for the semi-infinite lattice, the left-most site is occupied with probability

$$\pi_{1}=p_{\infty}=1-e^{-1}=0.6321205588\ldots$$

(16)

which substantially exceeds the probability (9) that a bulk site is occupied.

The probability $q_{L}$ that the second site, viz. the site adjacent to the left-most site is occupied, is dual to the probability $p_{L}$:

$$q_{L}=1-p_{L}=\sum_{k=0}^{L}\frac{(-1)^{k}}{k!}$$

(17)

for $L\geq 2$. In particular

$$\pi_{2}=q_{\infty}=e^{-1}=0.36787944117\ldots$$

(18)

The probability $r_{L}$ that the third site is occupied satisfies the recurrence

$$r_{L}=\frac{1}{L}+\frac{1}{L}\sum_{k=3}^{L-2}r_{k}+\frac{1}{L}\,p_{L-2}$$

(19)

which is established similarly to the recurrence (11). Using the generating function

$$R(x)=\sum_{L\geq 3}r_{L}\,x^{L}$$

(20)

we convert the recurrence (19) into a differential equation

$$\frac{dR(x)}{dx}=\frac{x}{1-x}\,R(x)+\frac{x^{2}}{1-x}+xP(x)$$

(21)

with $P(x)$ given by (14). Solving (21) subject to $R(0)=0$ we obtain

$$R(x)=\frac{1-x+x^{2}-\left(1+\frac{x^{2}}{2}\right)e^{-x}}{1-x}$$

(22)

from which

$$r_{L}=1-\frac{1}{2}\sum_{k=0}^{L-2}\frac{(-1)^{k}}{k!}-\sum_{k=0}^{L}\frac{(-1)^{k}}{k!}$$

(23)

In particular,

$$\pi_{3}=r_{\infty}=1-\tfrac{3}{2}e^{-1}=0.4481808382\ldots$$

(24)

A longer calculation gives

$$\pi_{4}=\tfrac{37}{4}e^{-1}-3=0.4028848308\ldots$$

(25)

These results hint that

$$\pi_{2n}=A_{n}e^{-1}-B_{n},\qquad\pi_{2n-1}=C_{n}-D_{n}e^{-1}$$

with some positive rational $A_{n},B_{n},C_{n},D_{n}$. Recalling that $\pi_{\infty}=\rho_{\text{jam}}$, one gets

$$\begin{split}&\lim_{n\to\infty}\left[1+2B_{n}-2A_{n}e^{-1}\right]=e^{-2}\\ &\lim_{n\to\infty}\left[1-2C_{n}+2D_{n}e^{-1}\right]=e^{-2}\end{split}$$

The total number of jammed states $J_{L}$ satisfies the recurrence

$$J_{L}=J_{L-2}+J_{L-3}$$

(26)

Indeed, the jammed states can be divided into the complementary sets with the occupied and empty leftmost site. If in a jammed state the leftmost site is occupied, the second site must be empty; the total number of such jammed states is equal to $J_{L-2}$. If the leftmost site in a jammed state is empty, the second site must be occupied and the third must be empty; the total number of such jammed states is equal to $J_{L-3}$.

The solution of the recurrence (26) subject to the boundary conditions $J_{0}=J_{1}=1$ and $J_{2}=2$ is encapsulated in the generating function

$$\sum_{L\geq 0}J_{L}x^{L}=\frac{1+x+x^{2}}{1-x^{2}-x^{3}}$$

(27)

In particular, the asymptotic behavior is

$$J_{L}\simeq\frac{(1+r)^{2}}{2r+3}\,r^{L}\quad\text{as}\quad L\to\infty$$

(28)

where $r$ is the real root of the cubic equation $r^{3}=r+1$:

	$\displaystyle r$	$\displaystyle=$	$\displaystyle\sqrt[3]{1+\sqrt[3]{1+\sqrt[3]{1+\ldots}}}$
		$\displaystyle=$	$\displaystyle\frac{\sqrt[3]{108+12\sqrt{69}}+\sqrt[3]{108-12\sqrt{69}}}{6}$
		$\displaystyle=$	$\displaystyle 1.3247179572\ldots$

known as the plastic number. The Padovan sequence (26) is a cousin of the Fibonacci sequence, and it also arises in several applications (see, e.g., Zagier ; JML ).

Consider an interval of length $L$. The RSA procedure brings the system into a jammed state. The total number $N$ of absorbed particles fluctuates from realization to realization. Here we compute the cumulant generating function that encodes all the cumulants of $N$. Thus we want to compute the average

$$F(\lambda,L)\equiv\langle e^{\lambda N}\rangle=\sum_{N}e^{\lambda N}P(N,L)$$

(29)

where $P(N,L)$ is the probability to have $N$ particles in a jammed state. The standard relation

$$\ln\langle e^{\lambda N}\rangle=\sum_{n\geq 1}\frac{\lambda^{n}}{n!}\,\langle N^{n}\rangle_{c}$$

(30)

then gives all the cumulants: the average $\langle N\rangle_{c}=\langle N\rangle$, the variance $\langle N^{2}\rangle_{c}=\langle N^{2}\rangle-\langle N\rangle^{2}$, etc.

The function $F(\lambda,L)\equiv\langle e^{\lambda N}\rangle$ grows exponentially with $L$. The cumulant generating function

$$U(\lambda)=\lim_{L\to\infty}L^{-1}\ln F(\lambda,L)$$

(31)

encapsulates all cumulants:

$$U(\lambda)=\sum_{n\geq 1}\frac{\lambda^{n}}{n!}\,U_{n},\qquad\langle N^{n}\rangle_{c}=LU_{n}$$

(32)

The cumulant generating function satisfies

$$F(\lambda,L)=\frac{e^{\lambda}}{L}\sum_{k=1}^{L}F(\lambda,k-2)F(\lambda,L-k-1)$$

(33)

To derive this recurrence, suppose $k$ is the first landing site. The intervals on the left and right of site $k$ are filled independently. (The one-dimensional nature of the problem is crucial for this property.) If $N_{-}$ and $N_{+}$ are the final occupation numbers of the left and right intervals, the total occupation number is $N=N_{-}+1+N_{+}$, so $e^{\lambda N}=e^{\lambda}e^{\lambda N_{-}}e^{\lambda N_{+}}$. Performing the averaging, summing over all $k$, and taking into account that site $k$ is chosen with probability $L^{-1}$, we obtain Eq. (33).

The recurrence (33) applies for all $L\geq 1$ if we set

$$F(\lambda,-1)=F(\lambda,0)=1$$

(34)

We now introduce the generating function

$$\Phi(\lambda,x)=\sum_{L\geq 0}F(\lambda,L)\,x^{L}$$

(35)

To recast (33) into an equation for the generating function we multiply (33) by $Lx^{L-1}$ and sum over all $L\geq 1$. The left-hand side turns into

$$\sum_{L\geq 1}LF(\lambda,L)\,x^{L-1}=\frac{\partial\Phi(\lambda,x)}{\partial x}$$

The right-hand side becomes

$$e^{\lambda}\sum_{L\geq 1}\sum_{k=1}^{L}F(\lambda,k-2)F(\lambda,L-k-1)\,x^{L-1}$$

and using the boundary conditions (34) we simplify the double sum to

$$x^{2}\sum_{m\geq-1}F(\lambda,m)\,x^{m}\sum_{n\geq-1}F(\lambda,n)\,x^{n}=\left[1+x\Phi(\lambda,x)\right]^{2}$$

and recast the recurrence (33) into a Riccati equation

$$\frac{d\Phi}{dx}=e^{\lambda}\left[1+x\Phi\right]^{2}$$

(36)

which turns out to be solvable. The solution of Eq. (36) subject to the initial condition $\Phi(\lambda,0)=1$ reads

$$\Phi(\lambda,x)=\frac{1+\Lambda\tanh(\Lambda x)}{1-x+(\Lambda^{-1}-\Lambda x)\tanh(\Lambda x)}$$

(37)

where $\Lambda\equiv e^{\lambda/2}$. The generating function $\Phi(\lambda,x)$ has a simple pole at $x=y=y(\lambda)$. Using (37) we find that the pole is implicitly determined by

$$\frac{\tanh(\Lambda y)}{\Lambda}=\frac{y-1}{1-\Lambda^{2}y}\,,\qquad\Lambda=e^{\lambda/2}$$

(38)

The cumulant generating function

$$U(\lambda)=-\ln y(\lambda)$$

(39)

is plotted in Fig. 1.

We already know the dominant exponential factor in $F(\lambda,L)\propto e^{LU(\lambda)}$. The exact result (37) allows one to find a pre-exponential factor. A straightforward calculation gives the residue of the simple pole of $\Phi(\lambda,x)$. We find $\Phi(\lambda,x)=(\Lambda y)^{-2}(y-x)^{-1}+O(1)$ when $x\to y$, from which we deduce

$$F(\lambda,L)\simeq(\Lambda^{2}y^{3})^{-1}\,e^{LU(\lambda)}$$

(40)

when $L\gg 1$.

The cumulant generating function is implicitly defined through (38)–(39). Expanding $U(\lambda)$, we recover the already known value of the first cumulant: $\langle N\rangle=L\rho_{\text{jam}}$. In principle, one can determine an arbitrary cumulant. We haven’t succeeded in finding an explicit general formula. Here we merely list a few cumulants that are extracted by expanding $U(\lambda)$ implicitly given by (38)–(39):

$$\begin{split}\frac{\langle N^{2}\rangle_{c}}{L}&=\frac{1}{e^{4}}\\ \frac{\langle N^{3}\rangle_{c}}{L}&=\frac{e^{4}-61}{16\,e^{6}}\\ \frac{\langle N^{4}\rangle_{c}}{L}&=\frac{43-e^{4}}{2\,e^{8}}\\ \frac{\langle N^{5}\rangle_{c}}{L}&=\frac{305e^{4}-10283-2e^{8}}{64\,e^{10}}\\ \frac{\langle N^{6}\rangle_{c}}{L}&=\frac{47871-1720e^{4}+17e^{8}}{32\,e^{12}}\\ \frac{\langle N^{7}\rangle_{c}}{L}&=\frac{359905e^{4}+17e^{12}-8540689-4697e^{8}}{512\,e^{14}}\\ \frac{\langle N^{8}\rangle_{c}}{L}&=\frac{6935126-335097e^{4}+5418e^{8}-31e^{12}}{32\,e^{16}}\end{split}$$

(41)

The so-called Mandel $Q$ parameter Mandel defined via

$$Q=\frac{\langle N^{2}\rangle_{c}}{\langle N\rangle}-1$$

(42)

is a basic measure characterizing the deviation from Poissonian statistics. The values $-1\leq Q<\infty$ are permissible, for the Poisson statistics $Q=0$ and the range $-1\leq Q<0$ is sub-Poissonian. The numerical value

$$Q=\frac{2}{e^{2}(e^{2}-1)}-1=-0.95763528097389\ldots$$

(43)

indicates that the statistics is strongly sub-Poissonian in the present case.

The ratios $\langle N^{n}\rangle_{c}/\langle N\rangle$ of cumulants to the average are known as Fano factors Fano . Here are a few Fano factors

$$\begin{split}\frac{\langle N^{3}\rangle_{c}}{\langle N\rangle}&=\frac{e^{4}-61}{8\,e^{4}(e^{2}-1)}=-0.0022940394167\ldots\\ \frac{\langle N^{4}\rangle_{c}}{\langle N\rangle}&=\frac{43-e^{4}}{e^{6}(e^{2}-1)}=-0.00449971626\ldots\\ \frac{\langle N^{5}\rangle_{c}}{\langle N\rangle}&=\frac{305e^{4}-10283-2e^{8}}{16\,e^{10}(e^{2}-1)}=0.00009049346\ldots\\ \frac{\langle N^{6}\rangle_{c}}{\langle N\rangle}&=\frac{47871-1720e^{4}+17e^{8}}{16\,e^{12}(e^{2}-1)}=0.0002787944\ldots\end{split}$$

For the Poisson distribution, all Fano factors are equal to unity. Thus Fano factors further illustrate a substantial deviation from the Poisson statistics.

Extremal jammed states are the states with the largest and the smallest number of occupied sites. Let us extract the probabilities of such states from the behavior of the cumulant generating function in the $\lambda\to\pm\infty$ limits. An asymptotic analysis of (38) yields (see also Fig. 1)

$$U(\lambda)=\begin{cases}\frac{1}{3}\lambda-\frac{1}{3}\ln 3&\lambda\to-\infty\\ \frac{1}{2}\lambda-\ln u&\lambda\to\infty\end{cases}$$

(44)

We have dropped the terms vanishing in the $\lambda\to\pm\infty$ limits and denoted by $u$ the positive root of the equation $u\,\tanh(u)=1$; numerically $u=1.19967864\ldots$. Combining (40) and (44) we get

$$F(\lambda,L)\sim\begin{cases}3^{-L/3}~{}e^{\lambda L/3}&\lambda\to-\infty\\ u^{-L}~{}e^{\lambda L/2}&\lambda\to\infty\end{cases}$$

(45)

To appreciate the first asymptotic in (45) we note that in the $\lambda\to-\infty$ limit the dominant contribution to the sum in Eq. (29) is provided by the jammed state with the smallest number of occupied sites. This jammed state

	$$\ldots\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\ldots$$		(46a)
has $N_{\text{min}}=\lfloor(L+2)/3\rfloor$ particles explaining the dominant $e^{\lambda L/3}$ factor. The $3^{-L/3}$ pre-factor gives the probability to end up in such jammed state.

In the $\lambda\to\infty$ limit the dominant contribution to the sum (29) is provided by the jammed state with the largest number of occupied sites:

$$\ldots\circ\,\bullet\,\circ\,\bullet\,\circ\,\bullet\,\circ\,\bullet\,\circ\,\bullet\,\circ\,\bullet\,\circ\,\bullet\,\circ\,\bullet\,\circ\,\bullet\,\circ\ldots$$

(46b)

We have $N_{\text{max}}=\lfloor(L+1)/2\rfloor$ and this explains the $e^{\lambda L/2}$ factor in (45). The $u^{-L}$ pre-factor is the probability to end up in the jammed state with the largest number of occupied sites. Thus the large deviation technique gives the asymptotic behaviors of the probability to reach the extremal jamming jammed states:

		$\displaystyle P(N_{\text{max}},L)\sim u^{-L}$		(47a)
		$\displaystyle P(N_{\text{min}},L)\sim 3^{-L/3}$		(47b)

The probabilities $P(N_{\text{max}},L)$ and $P(N_{\text{min}},L)$ can be also computed exactly as follows. Let us start with maximally dense jammed states (46b). The maximally dense state with $N_{\text{max}}=n$ particles may occur when $L=2n-1$ or $L=2n$. The probability $\mu_{n}=P(n,2n-1)$ can be recurrently determined from

$$\mu_{n}=\frac{1}{2n-1}\sum_{k=0}^{n-1}\mu_{k}\mu_{n-k-1}$$

(48)

Indeed, to ensure the maximal occupancy, the first landing must be to a site with an odd label, say $2k+1$. The interval of length $2k-1$ on the left (site $2k$ must remain unoccupied) and the interval of length $2(n-k-1)-1$ on the right are filled independently, and the requirement of the maximal occupancy gives the recurrence (48).

The boundary conditions read

$$\mu_{0}=\mu_{1}=1$$

(49)

Introducing the generating function

$$\mu(x)=\sum_{k\geq 0}\mu_{k}x^{k}$$

(50)

we convert the recurrence (48) into a Riccati equation

$$2\,\frac{d\mu}{dx}=\frac{\mu-1}{x}+\mu^{2}$$

(51)

which is solved to yield

$$\mu(x)=\frac{1}{1-\sqrt{x}\,\tanh\sqrt{x}}$$

(52)

This generating function has a simple pole $2/(u^{2}-x)$ at $x=u^{2}$, from which we extract the large $n$ asymptotic

$$\mu_{n}=P(n,2n-1)\simeq\frac{2}{u^{2n+2}}$$

(53)

The probability $\widehat{\mu}_{n}=P(n,2n)$ can be recurrently determined from

$$\widehat{\mu}_{n}=\frac{1}{n}\sum_{k=0}^{n-1}\widehat{\mu}_{k}\,\mu_{n-k-1}$$

(54)

which is derived similarly to (48). To appreciate the factor $n^{-1}$ in front of the sum we note that the probability of the first landing attempt is now $(2n)^{-1}$. The interval of even length may be either on the left or the right from the first landing site. This gives an extra factor of two.

The generating function

$$\widehat{\mu}(x)=\sum_{k\geq 0}\mu_{k}x^{k}$$

(55)

converts the recurrence (54) into a differential equation

$$\frac{d\widehat{\mu}}{dx}=\widehat{\mu}\,\mu$$

(56)

with known $\mu(x)$, see (52). Integrating (56) yields

$$\widehat{\mu}(x)=\frac{1}{\big{(}\cosh\sqrt{x}-\sqrt{x}\,\sinh\sqrt{x}\big{)}^{2}}$$

(57)

This generating function behaves as

$$\widehat{\mu}(x)=\frac{4(1-u^{-2})}{(u^{2}-x)^{2}}+\frac{2(1-u^{-2})}{u^{2}(u^{2}-x)}+O(1)$$

(58)

when $x\to u^{2}$, from which we extract the large $n$ asymptotic

$$\widehat{\mu}_{n}=P(n,2n)\simeq 2\left(1-u^{-2}\right)\frac{2n+3}{u^{2n}}$$

(59)

The asymptotic behaviors (53) and (59) are the more precise versions of the asymptotic (47a) following from the cumulant generating function.

Consider now the maximally sparse jammed states, say for intervals of length $L=3n$, so that $N_{\text{min}}=n$. The probabilities $\nu_{n}=P(n,3n)$ satisfy the recurrence

$$\nu_{n}=\frac{1}{3n}\sum_{k=0}^{n-1}\nu_{k}\nu_{n-k-1}$$

(60)

One can solve this recurrence using the generating function technique, from which one gets

	$$P(n,3n)=3^{-n}$$		(61a)
This simple solution can be verified by substitution of (61a) into the recurrence (60).

A similar analysis gives two other series of probabilities of maximally sparse jammed states:

		$\displaystyle P(n,3n-1)=\frac{3n+2}{5}\,\frac{1}{3^{n-1}}$		(61b)
		$\displaystyle P(n,3n-2)=\frac{63n^{2}+45n+8}{350}\,\frac{1}{3^{n-2}}$		(61c)

The exact results (61a)–(61c) are the more precise versions of the asymptotic (47b) following from the cumulant generating function.

Thus, in the realm of the RSA model, the defect-free arrays arise with exponentially small probabilities. The defect-free arrays (46a)–(46b) of Rb atoms excited to Rydberg states have been recently assembled experimentally by Bernien et al. Lukin17 .

The probability distribution $P(N,L)$ satisfies a recurrence relation

$\displaystyle LP(N,L)\!=\!\sum_{k=1}^{L}\sum_{i+j=N-1}P(i,k-2)P(j,L-k-1)$

(62)

which holds for all $L\geq 1$ if we set

$$P(N,L)=\delta_{N,0}\qquad\text{for all}\quad L\leq 0.$$

(63)

One can verify that $P(N,1)=\delta_{N,1}$ and $P(N,2)=\delta_{N,1}$. We now introduce the generating function

$$\mathcal{P}(x,y)=\sum_{N\geq 0}\sum_{L\geq 0}P(N,L)x^{N}y^{L}$$

(64)

This definition implies that

$$\sum_{N\geq 0}\sum_{L\geq 0}LP(N,L)x^{N}y^{L-1}=\frac{\partial\mathcal{P}}{\partial y}$$

(65)

Using (62)–(65) we obtain a differential equation for the generating function

$$\frac{\partial\mathcal{P}}{\partial y}=x[1+y\mathcal{P}]^{2}$$

(66)

which is mathematically identical to Eq. (36). The solution of (66) subject to the boundary condition $\mathcal{P}(x,0)=1$ is thus identical to (37) in different notation:

$$\mathcal{P}(x,y)=\frac{1+\sqrt{x}\tanh(y\sqrt{x})}{1-y+(\frac{1}{\sqrt{x}}-y\sqrt{x})\tanh(y\sqrt{x})}$$

(67)

The occupation number $N$ always lies within the bounds

$$\left\lfloor\frac{L+2}{3}\right\rfloor\leq N\leq\left\lfloor\frac{L+1}{2}\right\rfloor$$

(68)

For instance, $5\leq N\leq 7$ when $L=14$; expanding (67) one finds the corresponding probabilities

$$P(5,14)=\tfrac{17}{405},~{}~{}P(6,14)=\tfrac{1854632}{3274425},~{}~{}P(7,14)=\tfrac{1282348}{3274425}$$

The average number of occupied sites obeys

$$A_{L}=\frac{1}{L}\sum_{k=1}^{L}\big{(}A_{k-b-1}+1+A_{L-k-b}\big{)}$$

(69)

This recurrence is applicable for all $L\geq 1$ after setting $A_{j}=0$ for $j\leq 0$. Using the generating function (4) we convert (69) into a differential equation

$$\frac{dA}{dx}=\frac{2x^{b}}{1-x}\,A+(1-x)^{-2}$$

(70)

Solving (70) subject to $A(0)=0$ yields

$$A(x)=\frac{e^{-2L_{b}(x)}}{(1-x)^{2}}\int_{0}^{x}dy\,e^{2L_{b}(y)}$$

(71)

where we shortly write

$$L_{b}(x)=\sum_{j=1}^{b}\frac{x^{j}}{j}$$

(72)

Using (71) we extract the density of occupied sites

$$\rho_{b}=e^{-2L_{b}(1)}\int_{0}^{1}dy\,e^{2L_{b}(y)}$$

(73)

in the $L\to\infty$ limit.

For the minimal model we recover $\rho_{1}=(1-e^{-2})/2$. The next jammed density $\rho_{2}$ also admits an explicit expression through standard special functions:

$$\rho_{2}=\frac{\sqrt{\pi}}{2e^{4}}\left[\text{Erfi}(1)+\text{Erfi}(2)\right]=0.2745509877\ldots$$

(74)

where $\text{Erfi}(z)=-\sqrt{-1}\,\text{Erf}(\sqrt{-1}z)$ is an error function. The following jammed densities are

$$\rho_{3}\approx 0.200973,\quad\rho_{4}\approx 0.158455,\quad\rho_{5}\approx 0.130772$$

The maximally dense and maximally sparse jammed states provide obvious upper and lower bounds on the jammed density:

$$\frac{1}{2b+1}<\rho_{b}<\frac{1}{b+1}$$

(75)

These bounds hint that the limit

$$\lim_{b\to\infty}b\rho_{b}=C$$

(76)

exists and further suggest that $\frac{1}{2}<C<1$. One can compute $C$ using (73) and performing an asymptotic analysis. First we re-write (73) as

$$\rho_{b}=\int_{0}^{1}dy\,\exp\left[-2\sum_{j=1}^{b}\frac{1-y^{j}}{j}\right]$$

(77)

We then write $y=1-v/b$ and transform (77) into

$$b\rho_{b}=\int_{0}^{b}dv\,\exp\left[-2\sum_{j=1}^{b}\frac{1-\left(1-\frac{v}{b}\right)^{j}}{j}\right]$$

(78)

In the $b\to\infty$ limit, we can use the asymptotic relation $\left(1-\frac{v}{b}\right)^{j}\to e^{-vj/b}$ and replace summation over $j$ by integration over $u=vj/b$. This leads to

$\displaystyle C=\int_{0}^{\infty}dv\,\exp\left[-2\int_{0}^{v}du\,\frac{1-e^{-u}}{u}\right]$

known as the Rényi’s parking constant R58 .

The efficiency of the coverage is measured by the ratio of the jamming density $\rho_{b}$ to the maximal possible density: $\theta_{b}=(b+1)\rho_{b}$. This quantity monotonically decays as $b$ increases and it approaches $C=\theta_{\infty}$. Here are a few numerical values

$$\begin{split}\theta_{1}&=0.86466471676338\ldots\\ \theta_{2}&=0.82365296317734\ldots\\ \theta_{3}&=0.80389347991537\ldots\\ \theta_{4}&=0.79227591371305\ldots\\ \theta_{5}&=0.78463015586503\ldots\\ \theta_{\infty}&=0.74759791502876\ldots\end{split}$$

The function $F(\lambda,L)$ satisfies the recurrence

$$F(\lambda,L)=\frac{e^{\lambda}}{L}\sum_{k=1}^{L}F(\lambda,k-b-1)F(\lambda,L-k-b)$$

(79)

which we recast into a differential equation for the generating function $\Phi=\Phi(\lambda,x)$ defined by Eq. (35):

$$\frac{d\Phi}{dx}=e^{\lambda}\left[\frac{1-x^{b}}{1-x}+x^{b}\Phi\right]^{2}$$

(80)

This Riccati equation (80) appears exactly solvable only for the minimal $b=1$ model. It is still possible to determine the variance by analyzing $\partial_{\lambda}\Phi(\lambda,x)$ and $\partial_{\lambda}^{2}\Phi(\lambda,x)$ at $\lambda=0$. These functions satisfy equations which can be deduced from Eq. (80); by solving these equations one can determine the variance as we show in the appendix.

We also note that the probability distribution $P(N,L)$ satisfies a recurrence relation

$$LP(N,L)\!=\!\sum_{k=1}^{L}\sum_{i+j=N-1}P(i,k-b-1)P(j,L-k-b)$$

(81)

generalizing (62). Using the generating function (64) we recast the recurrence (81) into a Riccati equation

$$\frac{\partial\mathcal{P}}{\partial y}=x\left[\frac{1-y^{b}}{1-y}+y^{b}\mathcal{P}\right]^{2}$$

(82)

which is mathematically identical to (80) and appears unsolvable when $b=2,3,4,\ldots$.

The probabilities of the largest possible deviations, i.e. extremal jammed states, can be probed analytically using the same approach as in Sec. III.2. We start with maximally dense jammed states and choose for concreteness such states with $N_{\text{max}}=n$ in the chains of length $L=(b+1)n-b$. This

$$\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\bullet\,\circ\,\circ\,\,\bullet$$

(83)

is an example of such maximally dense jammed state with $N_{\text{max}}=6$ and $L=16$ for the model with the blockage radius is $b=2$. The probabilities $\mu_{n}=P(n,(b+1)n-b)$ satisfy the recurrence

$$\mu_{n}=\frac{1}{(b+1)n-b}\sum_{k=0}^{n-1}\mu_{k}\mu_{n-k-1}$$

(84)

and the boundary condition (49). The generating function (50) satisfies a Riccati equation

$$(b+1)\,\frac{d\mu}{dx}=b\,\frac{\mu-1}{x}+\mu^{2}$$

(85)

Solving this Riccati equation subject to $\mu(0)=1$ yields

$$\mu=-\frac{\sqrt{bx}\left[I_{\beta-1}(z)+I_{\beta-3}(z)\right]+(1+2b)I_{\beta-2}(z)}{2x\,I_{\beta-2}}$$

(86)

Here $I_{p}$ is the modified Bessel function and we have used shorthand notation $\beta=(1+b)^{-1}$ and $z=2\beta\sqrt{bx}$.

The denominator in (86) vanishes when $x\to y$, where $y=y(b)$ is found from

$$I_{\beta-2}\left(2\beta\sqrt{by}\right)=0,\qquad\beta=(1+b)^{-1}$$

(87)

This zero is simple and hence $\mu_{n}\sim[y(b)]^{-n}$. Re-writing in terms of the length of the chain, we arrive at

$$P(N_{\text{max}},L)\sim[u(b)]^{-L},\qquad u(b)=[y(b)]^{1/(b+1)}$$

(88)

The plot of $u(b)$ is shown in Fig. 2. Only positive integer values of the blockage radius, $b\in\mathbb{Z}_{+}$, matter. A few first numerical values are

$$\begin{split}u(1)&=1.19967864026\ldots\\ u(2)&=1.21712361233\ldots\\ u(3)&=1.20704834272\ldots\\ u(4)&=1.19255384190\ldots\\ u(5)&=1.17847101637\ldots\end{split}$$