Optimal random deposition of interacting particles

TL;DR

This paper derives exact solutions for the kinetics of interacting particle deposition in one dimension, revealing how to optimize packing density through interaction potentials, with implications for biological systems like DNA nucleosome arrangement.

Contribution

It provides the first exact analytical solution for one-dimensional interacting particle deposition with finite-range potentials, identifying a unique optimal potential for maximum density.

Findings

Exact time-dependent solutions for deposition processes

Existence of a unique optimal interaction potential

Optimal packing requires highly coordinated dynamics

Abstract

Irreversible random sequential deposition of interacting particles is widely used to model aggregation phenomena in physical, chemical, and biophysical systems. We show that in one dimension the exact time dependent solution of such processes can be found for arbitrary interaction potentials with finite range. The exact solution allows to rigorously prove characteristic features of the deposition kinetics, which have previously only been accessible by simulations. We show in particular that a unique interaction potential exists that leads to a maximally dense line coverage for a given interaction range. Remarkably, this distribution is singular and can only be expressed as a mathematical limit. The relevance of these results for models of nucleosome packing on DNA is discussed. The results highlight how the generation of an optimally dense packing requires a highly coordinated packing dynamics, which can be effectively tuned by the interaction potential even in the presence of intrinsic randomness.

Equations

$$\frac{\partial }{\partial t}p(x,t)$$

$$\varphi (t)=1-{\int }_{\sigma }^{\infty }\mathrm{d}x(x-\sigma )p(x,t),$$

$$\Omega (x,y)=k_{+}\omega (x)\omega (x-y),$$

$$\Omega (x,y)$$

$$\psi (x)$$

$$\frac{\partial }{\partial t}p$$

$$F(t)=\exp \left[-2{\int }_0^t\mathrm{d}s\frac{1-{\int }_{\Delta }^a\mathrm{d}u\rho (u)e^{-us}}{s}\right].$$

$$p_0(x,t)=t^{2}F(t)e^{-(x-2\sigma )t}.$$

$$\frac{\partial }{\partial t}{p}_1(x,t)$$

$$p_1(x,t)$$

$$p(x,t)=\{\begin{array}{cc}{p}_0(x,t),& x\ge 2a\\ & \\ p_1(x,t),& 2a-\Delta \le x<2a\\ ⋮& \\ p_j(x,t)& 2a-j\Delta \le x<2a-(j-1)\Delta \\ ⋮& \\ p_n(x,t)& 2\Delta \le x<2a-(n-1)\Delta \\ & \\ p_{n+1}(x,t)& \Delta <x<2\Delta \end{array}$$

$$p_j(x,t)$$

$$p_{n+1}(x,t)$$

$$p_{\mathrm{s}}(x)=2{\int }_0^{\infty }\mathrm{d}s{\int }_{x+\Delta }^{\infty }\mathrm{d}y\Omega (x,y)\sum _{i=0}^n{\Theta }_i(y)p_i(y,s).$$

$$\frac{\partial }{\partial t}{p}_{<}(x,t)$$

$$p_{<}(x,t)$$

$$h(x,t)$$

$$V(x)\approx (a-x)\kappa -\log \left[1+(a-x)\left(1-e^{-\kappa }\right)\right],$$

$$\Omega (x,y)=k_{+}\exp \left[-V(x)-V(y-x)\right],$$

$$p(x,t_{\mathrm{c}})=\{\begin{array}{cc}{p}_0(x,t_{\mathrm{c}}),& x\ge 2a\\ (1+\mathcal{L})h(x,t_{\mathrm{c}}),& 0\le x<2a\end{array}$$

$$V(x)=\{\begin{array}{cc}-\mu \log \left(\frac{x-\Delta }{a-\Delta }\right),& \Delta \le x\le a\\ 0,& x>a\end{array}$$

$${\tilde{p}}^{\ast }(z)=2{\int }_0^{\infty }\mathrm{d}tt\exp \left[-zt-2{\int }_0^t\mathrm{d}s\frac{1-e^{-as}}{s}\right],$$

$${\varphi }_{\mathrm{opt}}={\varphi }_{\mathrm{R}}+{\int }_{2\Delta -a}^a\mathrm{d}uz{\tilde{p}}^{\ast }(z).$$

$${\varphi }_{\mathrm{J}}-\varphi (t)\sim {\int }_{2\Delta }^{2a-(n-1)\Delta }\mathrm{d}x{e}^{-\psi (x)t}\sim t^{-\nu },$$

$$\psi (x)$$

$${\psi }_1(x)$$

$${\psi }_2(x)$$

$${\psi }_3(x)$$

$${\psi }_4(x)$$

$${\int }_0^t\mathrm{d}s\left(\frac{1-⟨e^{-us}⟩}{s}\right)$$

$$p_0(x,t)\simeq A_0{e}^{-(x-2\sigma )t},A_0=e^{2{\int }_0^{\infty }\mathrm{d}s⟨u{e}^{-us}⟩\log (s)}.$$

$$p_1(x,t)$$

Full text

The deposition of particles on a substrate is a ubiquitous phenomenon in science and engineering Elimelech (1995). From a theoretical perspective, the deposition dynamics is widely modelled in terms of a random sequential adsorption (RSA) process, which represents a paradigmatic adsorption mechanism Evans (1993). Since Renyi’s seminal work on the “car parking problem” (the RSA of equal line segments in 1$d$) Rényi (1958); Renyi (1963), RSA models and their variants have been successfully used to model polymer and colloid adsorption on surfaces Feder and Giaever (1980); Talbot et al. (2000); Cadilhe et al. (2007) and are also relevant in other contexts such as the compactification of granular matter Nowak et al. (1998); Tarjus and Viot (2004), genome sequencing Roach et al. (2000), and nucleosome packing on DNA Kornberg and Stryer (1988); Ranjith et al. (2007); Padinhateeri and Marko (2011); Möbius et al. (2013); Osberg et al. (2014). However, in almost all studies, particles have been assumed to interact solely through hard-core steric exclusion, despite the fact that realistic particles typically interact with a soft-core potential due to their internal structure Baule et al. (2018). In fact, including soft particle interactions significantly alters the filling behaviour and can account, e.g., for the observed rapid filling of nucleosomes on DNA Möbius et al. (2013); Osberg et al. (2014, 2015). In this work, we derive the exact time dependent solution of 1$d$ continuum RSA processes with arbitrary finite-range particle interactions using an iterative approach. This method also provides the exact general solution of the related RSA of polydisperse particles, which has been an open problem since the 1960s for the continuum 1$d$ case Ney (1962); Mullooly (1968); Krapivsky (1992); Brilliantov et al. (1996, 1997); Burridge and Mao (2004), notwithstanding specific scaling solutions that have been found Krapivsky (1992); Brilliantov et al. (1996, 1997); Hassan (1997). The exact solution allows in particular to address the fundamental question how the packing density of the deposition process can be optimized by tuning the particle interactions or size distribution. Remarkably, the solution reveals that a unique interaction potential/size distribution exists that leads to a maximally dense coverage of the line, which is approached as $\sim t^{-\nu}$, where $\nu\to 0^{+}$ is infinitesimally small.

In RSA, a particle’s position is selected with uniform probability over the domain and it is then placed sequentially if there is no overlap with any previously placed particles. Particles are not able to move or reorient once being placed. In order to parametrize the RSA dynamics of equal particles by a rate equation in $1d$, we define an interval $x$ as the distance between the centers of two nearest-neighbor particles on the line of length $L$ and introduce $N(x,t)\mathrm{d}x$ as the number of intervals with size $\in[x,x+\mathrm{d}x]$ at time $t$. The definition of $x$ implies that $\int_{0}^{L}\mathrm{d}x\,x\,N(x,t)=L$ is a conserved quantity of the dynamics for all $t$. The interval distribution $p(x,t)$ is defined as $p(x,t)=\frac{\lambda}{L}N(x,t)$, where $\lambda$ is a characteristic length scale associated with the particles 111For hard particles $\lambda$ would be identified as the size of the particles. For interacting (soft) particles, $\lambda$ could be identified as the interaction range, effective size, or hard core of the particles.. Scaling length as $x\to x/\lambda$ and considering the limit $L\to\infty$, the time evolution of $p(x,t)$ is exactly described by the master equation

[TABLE]

where the first term on the rhs describes the destruction of intervals of length $x$ and the second term the creation. $\Omega(x,y)$ is the probability per unit time that a particle is placed inside an interval of length $y$, thus creating an interval of length $x$. The factor 2 stems from the fact that in 1$d$ there are always two ways of doing this for a given $y$ ($x$ is either the distance of the newly inserted particle with the existing left or right particle). The destruction term is accordingly governed by the function $\psi(x)=\int_{0}^{x}\mathrm{d}u\,\Omega(u,x)$. The interval distribution satisfies $\int_{0}^{\infty}\mathrm{d}x\,x\,p(x,t)=1$ for all $t$. Moreover, since one particle is associated with every interval $x$, the integral $\int_{0}^{\infty}\mathrm{d}x\,p(x,t)=n(t)$ equals the number density of particles. The requirement of an initially empty line thus leads to the initial condition $n(0)=0$ or $p(x,0)=0$. From Eq. (1) this implies that $\lim_{t\to 0}\frac{\partial}{\partial t}p(x,t)=0$, but one can show that $\lim_{t\to 0}\frac{\partial^{2}}{\partial t^{2}}p(x,t)>0$ such that $n(t)$ monotonically increases with time. The key quantity of interest is the line coverage (packing density)

[TABLE]

where $\sigma$ denotes the effective size (diameter) of a particle (see below). Note that for soft particles $n(t)$ and $\phi(t)$ have very different long-time behaviours: while the jamming density $\phi_{\rm J}=\lim_{t\to\infty}\phi(t)\leq 1$, $n(t)$ can diverge as $t\to\infty$, since particles can be absorbed inside existing particles. By contrast, for hard particles they are directly related as $\phi(t)=\sigma n(t)$ Talbot et al. (2000).

The simplest example of a deposition process governed by Eq. (1) is Renyi’s seminal car parking problem, where particles only interact by steric repulsion Rényi (1958); Renyi (1963); Mackenzie (1962); Widom (1966). In this case, an exact solution for $p(x,t)$ is known and yields as hallmark features $\phi_{\rm J}=0.7475...\equiv\phi_{\rm R}$ with the algebraic asymptotic approach $\phi_{\rm J}-\phi(t)\sim t^{-1}$. Other models of the form Eq. (1) have been solved for specific $\Omega$ that lead to scale invariant solutions $p$, e.g, RSA Krapivsky (1992); Brilliantov et al. (1996, 1997) and fragmentation processes Ziff and McGrady (1985); Cheng and Redner (1988, 1990); Williams (1990). In the following, an exact analytical solution is derived much more generally without the requirement of scale invariance.

For interacting particles we assume the form

[TABLE]

where the adsorption rate $k_{+}$ sets the time scale (set to unity) and $\omega(x)$ describes the modification in the rate due to particle overlap. The particle interactions are constrained as follows: (i) $\omega(x)=1$ for $x\geq a$, i.e., $a$ is the finite range of the interaction; and (ii) $\omega(x)=0$ for $x\leq\Delta$, i.e., $\Delta$ is the hard-core exclusion volume of a particle. Clearly, (i,ii) are satisfied for almost all realistic particle models at least to a good approximation. The hard particle case of Renyi is obtained as $\omega(x)=\Theta(x-a)$ with $\Delta=a$, where $\Theta(x)$ denotes the Heaviside step function. For soft particles $\omega(x)$ interpolates between 0 and 1 (see Fig. 2a).

The properties (i,ii) of $\omega(x)$ constrain the form of $\psi$ and $\Omega(x,y)$ as follows (Appendix Sec. A):

[TABLE]

where $\sigma=a-\int_{\Delta}^{a}\mathrm{d}u\,\omega(u)$ can be interpreted as the effective size of a particle.

With Eqs. (4–6) the solution of Eq. (1) can be found as follows. For $x\geq 2a$, Eq. (1) simplifies to

[TABLE]

Crucially, Eq. (7) admits an exact closed form solution. We make the ansatz $p(x,t)=A(t)e^{-(x-2\sigma)t}$, which, upon substitution into Eq. (7), yields an equation for $A$: $\dot{A}/A=2\int_{\Delta}^{\infty}\mathrm{d}u\,\omega(u)e^{-ut}$. In order to make the divergence in the integral at $t=0$ explicit, we perform a partial integration: $\int_{\Delta}^{\infty}\mathrm{d}u\,\omega(u)e^{-ut}=t^{-1}\int_{\Delta}^{a}\mathrm{d}u\rho(u)e^{-ut}$, where the function $\rho$ is defined as $\rho(x)=\omega^{\prime}(x)$, which has support on $[\Delta,a]$ and satisfies $\int_{\Delta}^{a}\mathrm{d}u\,\rho(u)=1$. The boundary conditions further suggest the form $A(t)=t^{2}F(t)$ with $F(0)=1$ and we obtain the regularized form

[TABLE]

Let us denote the resulting solution of Eq. (7) as $p_{0}$, valid for $x\geq 2a$:

[TABLE]

Eq. (9) is valid for arbitrary $\Delta\geq 0$. We now distinguish the two cases $\Delta>0$ and $\Delta=0$.

The key observation is that for a finite excluded volume $\Delta>0$, we directly obtain a solution in the regime $2a-\Delta\leq x\leq 2a$ (denoted by $p_{1}$) even though $\psi$ is non-linear: property (ii) of $\omega(x)$ enforces the constraint $\Theta(y-x-\Delta)$ in the integral in Eq. (1), which thus only integrates over $p_{0}$, i.e., the master equation for $2a-\Delta\leq x\leq 2a$ is

[TABLE]

This simple first-order ODE with inhomogeneity can be directly integrated:

[TABLE]

Likewise, the solution in the range $2a-2\Delta\leq x\leq 2a-\Delta$ can be obtained by integration over $p_{0}$ and $p_{1}$, and so on. Overall, we thus decompose the interval distribution as

[TABLE]

Due to the excluded volume, $p(x,t)=0$ for $x\leq\Delta$. We introduce the shorthand notation $\Theta_{j}(x)$ to separate the different $x$ ranges in Eq. (12) and write $p(x,t)=\sum_{j=0}^{n+1}\Theta_{j}(x)p_{j}(x,t)$. The solutions $p_{1},...,p_{n+1}$ are thus obtained by solving the ODE Eq. (1) for each range leading to

[TABLE]

and

[TABLE]

since $\psi(x)=0$ for $x\leq 2\Delta$. Eqs. (9–Optimal random deposition of interacting particles) represent the exact analytical solution of Eq. (1) for finite range interactions and $\Delta>0$.

From Eq. (8) we obtain the asymptotic behaviour $p_{0}(x,t)\sim A_{0}e^{-(x-2\sigma)t}$ and in general $p_{j}(x,t)\sim A_{j}(x)e^{-\psi(x)t}$ for $j=1,...,n$, where the $A_{j}$ can be calculated iteratively (Appendix Sec. B). We see that the distributions $p_{1},...,p_{n}$ all decay exponentially governed by $\psi(x)$. The stationary limit is thus entirely determined by $p_{n+1}$: $\lim_{t\to\infty}p(x,t)=\lim_{t\to\infty}p_{n+1}(x,t)=p_{\rm s}(x)$, where from Eq. (13)

[TABLE]

With Eq. (15) the jamming density can be directly evaluated: $\phi_{\rm J}=1-\int_{\sigma}^{2\Delta}\mathrm{d}x(x-\sigma)p_{\rm s}(x)$.

When $\Delta\to 0$, we can use a limiting procedure in Eqs. (13,Optimal random deposition of interacting particles) or use Eq. (9) in Eq. (1) and separate the integration region. In both cases, the result is

[TABLE]

where $p_{<}$ is the interval distribution for the whole range $0\leq x\leq 2a$. Using the Duhamel principle, Eq. (16) can be solved by iteration leading to a series solution. We define the compact linear operator $\mathcal{L}f(x,t)=2\int_{0}^{t}\mathrm{d}s\,e^{-\psi(x)(t-s)}\int_{x}^{2a}\mathrm{d}y\,\Omega(x,y)f(y,t)$ and obtain the formal solution

[TABLE]

Convergence of this series needs to be established for a given $x,t$ range and $\Omega$.

As an example, where these results can provide novel analytical insight, we consider a model of nucleosome packing on DNA. In Möbius et al. (2013); Osberg et al. (2014) genome packaging in eukaryotic cells has been modelled assuming an effective ‘softness’ of the nucleosomes due to a multitude of internal states with different footprints on DNA Möbius et al. (2013); Osberg et al. (2014). The effective interaction potential has been approximated from in vivo data as Möbius et al. (2013)

[TABLE]

where $a=167$ base pairs (bp) is the maximal footprint size and $\kappa=0.15$ the stiffness per bp. A numerical investigation of the equilibrium kinetics of Eq. (19) exhibits a universal interval distribution independent of the adsorption rate $k_{+}$ at a specific time $t_{\rm c}$. The time $t_{\rm c}$ is defined by the condition $\phi(t_{\rm c})=\phi_{\rm R}$ and denotes the onset of ‘cramming’, where nucleosomes are increasingly squeezed into gaps that are smaller than the maximal footprint size Möbius et al. (2013); Osberg et al. (2014). The dynamics in the cramming regime can be described as follows. In nucleosome packing both adsorption (with rate $k_{+}$) and desorption (with rate $k_{-}$) occur whereby $k_{+}\gg k_{-}$. Thus in a time regime $\ll 1/k_{-}$ the equilibrium dynamics is well described by an irreversible RSA process like Eq. (1) Epstein (1979a, b). Here, the effect of the interaction potential can be captured by the Boltzmann factor Möbius et al. (2013); Osberg et al. (2014)

[TABLE]

setting $k_{B}T=1$ and the normalization constant is included in $k_{+}$. Since $t_{\rm c}$ has been observed to be $\ll 1/k_{-}$, the RSA process applies in the cramming regime and explains immediately the apparent universality of the interval distribution: different $k_{+}$ lead to different $t_{\rm c}$ values, but the resulting $p(x,t_{\rm c})$ are all identical, since they are determined at the same $\phi$ value. In fact, this observation shows that the universality holds not only at cramming onset, but for any fixed $\phi$ value in the regime $\ll 1/k_{-}$. It also highlights that the curves $\phi(t)$ for different $k_{+}$ are all scaled versions of each other up to times $t\approx 1/k_{-}$, which is indeed suggested in the numerical results of Osberg et al. (2014). Remarkably, with Eqs. (9,17,Optimal random deposition of interacting particles) we obtain

[TABLE]

which shows perfect agreement with the numerically obtained interval distribution at $t_{\rm c}$ Osberg et al. (2014) despite the strong nonlinear character of the interaction potential Eq. (19).

We now want to understand how dense packings on the line can be generated by tuning the interaction potential $V$ assuming $\omega(x)=e^{-V(x)}$ as in Eq. (20). For simplicity, we restrict the discussion to purely repulsive interactions such that $\omega(x)$ is monotonically increasing with $0\leq\omega(x)\leq 1$ for $x\in[\Delta,a]$. We then make two key observations: (a) For $\sigma\geq 2\Delta$ the effective size is larger than the minimal separation of two particles, thus the line will eventually be fully covered by particles and $\phi_{\rm J}=1$. We thus assume the more interesting case $2\Delta>\sigma$ in the following. (b) For $a>\Delta$ a potential leading to a maximum in $\phi_{\rm J}$ must exist. We can conclude this surprising fact from the two limiting forms of $\omega$ that satisfy the properties (i,ii) and the repulsive nature (see Fig. 2a). One limit is $\omega(x)=\Theta(x-a)$ (leading to $\sigma=a$) and the other limit is $\omega(x)=\Theta(x-\Delta)$ (leading to $\sigma=\Delta$). However, in both cases we recover $\phi_{\rm J}=\phi_{\rm R}$, since $\phi_{\rm J}$ is invariant with respect to a single size scale $a$ or $\Delta$ on an infinite line. When $a>\Delta$, we have $\phi_{\rm J}>\phi_{\rm R}$, so a potential leading to a maximum must exist for a given $a$. What is the form of this optimal potential?

In order to elucidate this matter, we consider a one-parameter family of potentials that can interpolate between the two limiting forms of $\omega(x)$ (see Fig. 2a)

[TABLE]

where $\mu>0$. The resulting $\phi_{\rm J}$ increases monotonically for larger $\mu$ (Fig. 2c). Surprisingly, this suggests that the maximally dense packing is reached when $\omega(x)$ becomes infinitesimally close to the step function limit $\Theta(x-a)$. One can show (Appendix Sec. C) that the interval distribution then converges to the stationary limit $p_{\rm s}(x)=\Theta(x-a)\Theta(2\Delta-x)\tilde{p}^{*}(x-a)$, where

[TABLE]

and the corresponding maximal line coverage is

[TABLE]

As shown in Fig. 2b, $p_{\rm s}(x)$ diverges at $x=\sigma=a$ indicating a large number of particle configurations with no empty space between neighbors. The optimal density $\phi^{\rm opt}$ increases monotonically as a function of $a$ (see inset of Fig. 2c). Note that the limit $\mu\to\infty$ is singular, recovering the Renyi density $\phi_{\rm R}$ instead (Appendix Sec. C). Remarkably, Eqs. (23,24) remain valid for any potential for which $\omega(x)$ approaches $\Theta(x-a)$ infinitesimally closely, highlighting that these results are independent of the specific form Eq. (22). This non-trivial property can be understood intuitively. In the $\mu\gg 1$ regime the system behaves initially like the Renyi car parking of hard particles: since $\omega(x)$ is infinitesimally close to $\Theta(x-a)$ all trial configurations with particle overlap are practically rejected until no gaps $x\geq 2a$ remain. Subsequently, the steep decay of $\omega(x)$ prevents that intervals $x<2a$ are filled with larger than necessary overlaps. The optimal potential thus induces a perfectly hierarchical filling, where all intervals $2\Delta<x<2a$ are eventually filled, but extremely slowly. In fact, with the previous results and the asymptotic properties of $\psi$ we find that for $\mu\gg 1$ (Appendix Sec. D)

[TABLE]

where $\nu=\frac{1}{2\mu+1}$. Nevertheless, Eq. (24) can be verified in a simulation by running the Renyi car parking RSA with length scale $a$ and then considering all intervals $x>2\Delta$ as filled when evaluating $\phi_{\rm J}$. This indeed yields perfect agreement with the theory (see inset of Fig. 2c).

The iterative solution method can be applied to related interacting particle models such as the continuum RSA of polydisperse particles in 1$d$, whose general time-dependent solution has been a long-standing open problem Ney (1962); Mullooly (1968); Krapivsky (1992); Brilliantov et al. (1996, 1997); Burridge and Mao (2004). Crucially, this problem can be solved in full generality akin to the interacting particle case for an arbitrary size distribution $\rho(\sigma)$ with support $\sigma\in[\sigma_{1},\sigma_{2}]$ (Appendix Sec. E). The solution shows in particular that the exponent of the asymptotic approach is determined by the leading term in the asymptotic expansion of $\rho(\sigma)$ as $\sigma\to\sigma_{1}$, which makes sense intuitively because the large time asymptotics is determined by the filling of the smallest available intervals. The widely used assumption that polydispersity generally reduces the exponent by an additional degree of freedom from $\nu=1/d$ to $\nu=1/(d+1)$ Tarjus and Talbot (1991); Evans (1993); Adamczyk et al. (1997); Talbot et al. (2000), where $d$ is the spatial dimension, is thus not correct for $d=1$. A similar result should hold for the RSA of polydisperse spheres in higher dimensions, where the quantitative dependence on $\rho$ has first been observed empirically Meakin and Jullien (1992). Remarkably, the size distribution $\rho(x)$ that leads to a maximal line coverage yields the same interval distribution and optimal density as in the interacting particle case, given by Eqs. (23,24), where now $\sigma_{1}=2\Delta-a$ and $\sigma_{2}=a$ (Appendix Sec. V).

The RSA of interacting particles coresponds to a cooperative RSA problem, where adsorption rates depend on the local environment of the particle with a given range Evans (1993). Cooperative RSA problems have long been studied on discrete lattices, for which exact results are available for $M$-mers (particles occupying $M$ sites) with cooperativity of range-1 González et al. (1974); Boucher (1973); Epstein (1979a), range-$M$ Wolf et al. (1984); Mellein (1986), and arbitrary finite range-$R$ Evans (1990). From such discrete models, the continuum equivalent is obtained in the limit of $M\to\infty$. However, performing such a limit in the general case of $M$-mers with general range-$R$ cooperativity, which could be mapped onto the continuum RSA of interacting particles, is non-trivial and left for future work. Further applications of the exact results in the context of nucleosome packing models would be highly interesting. The empirically obtained potential Eq. (19), e.g., might represent an optimal trade-off between achieving a high coverage and a fast filling dynamics on time scales relevant for biological function. To this end not only $\phi_{\rm J}$ but the dynamics of $\phi(t)$ needs to be explored in the space of possible potentials. In Gutiérrez et al. (2015), cold atoms that are excited to high-lying Rydberg states are shown to fill up excitation levels akin to a classical deposition process, where atoms interact through a highly nonlinear potential with excluded volume Lesanovsky and Garrahan (2013). The analytical solution obtained here could provide important deeper insight into such cold atom excitations, which would be readily realizable in experiments.

Appendix A Constraints on the master equation for finite-range interactions

Including explicitly the constraints (i,ii) on $\omega(x)$, we can express $\psi(x)=\int_{0}^{x}\mathrm{d}y\,\omega(x-y)\omega(y)$ as

[TABLE]

Factorization yields four terms

[TABLE]

We see immediately that $\psi(x)=0$ for $x\leq 2\Delta$, i.e., Eq. (5) holds. On the other hand, for $x\geq 2a$, the first term in Eq. (27) vanishes and the second and third terms both become $\int_{\Delta}^{a}\mathrm{d}y\,\omega(y)$, which leads to the linear form Eq. (6).

Appendix B Asymptotic time regime

We determine successively expressions for $p_{0},p_{1},...$ that are valid as $t\to\infty$. First, using partial integration in Eq. (8) we obtain

[TABLE]

where we use the shorthand notation $\langle...\rangle=\int_{\Delta}^{a}\mathrm{d}u\,...\rho(u)$. Asymptotically, $p_{0}$ of Eq. (9) becomes for large $t$

[TABLE]

In order to evaluate $p_{1}$, we note that Eq. (13) for $j=1$ can be simplified since the integrations only act on $p_{0}$ and we can substitute Eq. (7)

[TABLE]

If we introduce a time $t_{\rm a}$ such that the asymptotic behaviour Eq. (29) is valid for $t>t_{\rm a}$, we obtain from Eq. (30) by separating the integral

[TABLE]

We thus obtain the asymptotic behaviour for $t\to\infty$

[TABLE]

From Eq. (13) we obtain likewise for $p_{2}$

[TABLE]

The only problematic term is here the integral

[TABLE]

whose precise asymptotics for $t\to\infty$ can not be determined in general, since $\Omega(x,y)$ vanishes exponentially fast at $y=x+\Delta$, the minimum of $\psi(y)$. We assume that this term decays faster than $e^{-\psi(x)t}$, since $\psi(x+\Delta)>\psi(x)$ in general, but this needs to be argued for each potential $V$. With this assumption we obtain for $t>t_{\rm a}$

[TABLE]

The calculation can be continued for general $p_{j}$, $j=3,...,n$, which yields

[TABLE]

where $A_{j}$ depends on $A_{0},A_{1},...,A_{j-1}$. Even though the derivation of relies on assumptions of the decay of integrals of the type of Eq. (34), Eq. (37) is also motivated heuristically from the filling dynamics of Eq. (1): larger intervals are more easily filled than smaller ones, such that the $p_{j}$ need to decay hierarchically for $t\to\infty$. Since for finite $\Delta$ the creation term in Eq. (1) depends on $p_{0},p_{1},...,p_{j-1}$, it thus vanishes before $p_{j}$, leading to the asymptotic behaviour Eq. (37).

Appendix C Analysis of the $\mu\gg 1$ regime

Since we are interested in the jamming density $\phi_{\rm J}=\lim_{t\to\infty}\phi(t)$, which is determined from Eq. (2), we need to clarify the role of the two limits $t\to\infty$ and $\mu\to\infty$. First we note that for the potential Eq. (22) the following limits apply:

[TABLE]

In both limits the resulting jamming density is then just the Renyi density $\phi_{\rm R}$, since the particles are effectively hard with only a single length scale. For $\mu\to 0$ this follows formally from Eqs. (2,12,Optimal random deposition of interacting particles) as

[TABLE]

where the last step follows since for $\omega(x)=\Theta(x-\Delta)$ all distributions $p_{1},...,p_{n}$ reduce to $p_{0}$. Using $\lim_{\mu\to 0}p_{0}(x,t)=t^{2}\exp\left[-2\int_{0}^{t}\mathrm{d}s\frac{1-e^{-\Delta\,s}}{s}\right]e^{-(x-2\Delta)t}$, see Eq. (9), we obtain

[TABLE]

For the limit $\mu\to\infty$, we need to be more careful with the sequence of limits to recover $\phi_{\rm R}$. In particular, if we take $\mu\to\infty$ before $t\to\infty$, the distributions $p_{1},...,p_{n}$ no longer decay to zero for $t\to\infty$, because $\psi(x)=0$ for $x\leq 2a$. In fact, Eqs. (13) then take the form of Eq. (Optimal random deposition of interacting particles), i.e., $p_{1},...,p_{n}$ recover $p_{n+1}$ whose range of validity is now extended to the interval $[\Delta,2a]$. We thus obtain

[TABLE]

From then on the same calculation as above in Eq. (41) applies, where $\Delta$ is now replaced by $a$ leading also to $\phi_{\rm R}$. What happens now, when $\mu\gg 1$ but without taking the actual limit $\mu\to\infty$? For any large but finite $\mu$, the distributions $p_{1},...,p_{n}$ will eventually decay to zero for $t\to\infty$, since $\psi(x)$ is vanishingly small but non-zero for $x\leq 2a$. Apart from this crucial difference, the other limits Eqs. (38) hold effectively. As a result, Eq. (2) becomes

[TABLE]

Here, $\tilde{p}^{*}(z)$ is given in Eq. (23). We see that the stationary distribution is given by

[TABLE]

Noting that $\Theta(a-z)\tilde{p}^{*}(z)$ is the distribution of empty intervals in the Renyi car parking problem (see Eq. (C) with $\Delta=a$), the jamming density can be further expressed as

[TABLE]

The optimal jamming density $\phi^{\rm opt}$ is thus achieved for arbitrarily large but finite $\mu\gg 1$. As discussed above, the limit $\mu\to\infty$ is singular, recovering instead the Renyi density $\phi_{\rm R}$.

Appendix D Asymptotic approach to the jamming density

The asymptotic approach to $\phi_{\rm J}$ can be calculated with the large $t$ property Eq. (37). With Eq. (2), we obtain

[TABLE]

for large $t$, using Eq. (37) and the fact that $p_{n+1}$ converges to the stationary interval distribution, which determines $\phi_{\rm J}$. It remains to determine the behaviour of $\psi(x)$ in the vicinity of the minimum of the integration range at $x=2\Delta$. From Eq. (27) we obtain as $x\to 2\Delta$

[TABLE]

for the potential Eq. (22). Since there are no other singularities in the integrand of Eq. (46), we obtain the asymptotic scaling

[TABLE]

Appendix E Exact solution of the polydisperse RSA problem

We now consider a different RSA problem in 1$d$, where the particles are hard, but their sizes (diameters) $\sigma$ are drawn from an independent distribution $\rho(\sigma)$ at each adsorption attempt. We assume the PDF $\rho$ normalized with finite support $\sigma\in[\sigma_{1},\sigma_{2}]$. In this case, the time evolution of the interval distribution $p(x,t)$ is not described by Eq. (1), since the particles bounding every interval are of different size, which is not taken into account in Eq. (1). Formally, one would need to extend $p(x,t)$ to $p(x,t;\sigma_{\rm L},\sigma_{\rm R})$, where $\sigma_{\rm L},\sigma_{\rm R}$ denotes the sizes of the left and right particles separated by $x$. The time evolution of $p(x,t;\sigma_{\rm L},\sigma_{\rm R})$ follows a more complicated master equation, which can be mapped to the one of the ‘Paris car parking problem’ [35]. Nevertheless, we can describe the polydisperse RSA problem by an equation of the type of Eq. (1), if we focus on the distribution of empty gaps, $\tilde{p}(z,t)$, instead of the interval distribution $p(x,t)$ ($x$ is the distance between the centers of two particles, while $z$ is the uncovered part of $x$). The master equation (1) is then

[TABLE]

where the brackets $\langle...\rangle$ denote an expected value with respect to the size distribution $\rho$. The destruction term takes into account the possible ways of inserting a particle with size $\sigma$ sampled from $\rho$ in an empty gap of length $z$. Likewise the creation term is averaged over all possible sizes. As before, the rate constant is absorbed in the time scale and set to unity. We rewrite

[TABLE]

and, by comparing with Eqs. (1,4), identify

[TABLE]

where due to normalization $\omega(z)=1$ for $z\geq\sigma_{2}$. In analogy to Eqs. (5,6), $\tilde{\psi}(z)$ satisfies

[TABLE]

where $\overline{\sigma}=\left<\sigma\right>$. As a consequence, the same solution method as in the interacting particle case can be applied. The results equivalent to Eqs. (9–Optimal random deposition of interacting particles) are

[TABLE]

where

[TABLE]

and the line coverage is given by

[TABLE]

The $t\to\infty$ behaviour of $\tilde{p}_{1},...,\tilde{p}_{n}$ can be derived in analogy to the discussion in Sec. B showing again an exponential decay governed by $e^{-\psi(z)t}$.

Further below, we need the following expression for the rate of change $\dot{\phi}(t)$. With Eqs. (49,59) we obtain

[TABLE]

where we define $\Lambda(z)=\left<\sigma(z-\sigma)\Theta(z-\sigma)\right>$ in the last step.

E.1 Power-law size distribution

In order to find an optimal size distribution that leads to a maximally dense polydisperse packing, we follow the same approach as in the interacting particle case and consider a power-law size distribution

[TABLE]

such that the corresponding $\omega(z)$ is identical to the one resulting from the potential Eq. (22). However, it is important to note that the resulting dynamics of $\tilde{p}$ is different from that of the interacting particle case for finite $t$. The reason is that $\omega(x)=0$ for $x\leq\Delta$, but $\psi(x)=0$ for $x\leq 2\Delta$. In the polydisperse case, we have $\omega(z)=0$ and $\psi(z)=0$ for $z\leq\sigma_{1}$. However, the stationary distribution is identical in both cases for $\mu\gg 1$. With Eq. (58) we obtain

[TABLE]

For $\mu\gg 1$ we have $\omega(z)\to\Theta(z-\sigma_{2})$, $\sigma\to\sigma_{2}$, $\langle e^{-\sigma t}\rangle\to e^{-\sigma_{2}t}$ such that

[TABLE]

where with Eq. (56)

[TABLE]

Therefore, the stationary distribution of empty gaps is given by

[TABLE]

using Eq. (23) with $\sigma_{2}=a$. With Eq. (59) we obtain the jamming density

[TABLE]

which agrees with Eq. (24) provided we also set $\sigma_{1}=2\Delta-a$.

E.2 The case $n=1$

For the special case $\sigma_{1}<\sigma_{2}\leq 2\sigma_{1}$ a result for the jamming density $\phi_{\rm J}$ has been obtained in Ref. Burridge and Mao (2004) using a recursive approach. We reproduce this result in the following. With Eq. (60) we have the line coverage as a function of time

[TABLE]

calculating the $z$-integral in the first term. The distribution $\tilde{p}_{1}$ is given by Eq. (57). Using Eq. (52) we obtain

[TABLE]

Substituting $\tilde{p}_{1}$ in Eq. (67) and taking the $t\to\infty$ limit yields

[TABLE]

after a partial integration, which is the expression obtained in Burridge and Mao (2004).

E.3 Asymptotic scaling

Since the $\tilde{p}_{j}$ for $j=0,...,n-1$ all decay exponentially as $t\to\infty$, the asymptotic approach to $\phi_{\rm J}$ is governed by $\tilde{p}_{n}(z,t)\simeq B_{n}(z)e^{-\psi(z)t}$ in analogy to Sec. B. The scaling of $\phi_{\rm J}-\phi(t)$ could thus be calculated as in Sec. D, but here we use instead Eq. (60) to obtain

[TABLE]

As $z\to\sigma_{1}$, the minimum of $\psi(z)$, we have $\Lambda(z),\psi(z)\to 0$, while $B_{n}(z)$ remains finite (e.g., for $n=1$, $B_{n}=\tilde{p}_{1}(z,t_{\rm c})e^{\psi(z)t_{\rm c}}$). The asymptotic scaling thus depends on the behaviour of $\Lambda(z)/\psi(z)$ as $z\to\sigma_{1}$. We obtain

[TABLE]

and thus

[TABLE]

Therefore

[TABLE]

assuming the asymptotic form

[TABLE]

Since $\psi(\sigma_{1})=\psi^{\prime}(\sigma_{1})=0$ and $\psi^{\prime\prime}(z)=\rho(z)$, the exponent $\alpha$ is determined by the leading term in the asymptotic expansion of $\rho$ as $z\to\sigma_{1}$. For the general power-law size distribution Eq. (61) we obtain

[TABLE]