Exact hopping and collision times for two hard discs in a box

TL;DR

This paper derives exact analytical formulas for the mean times between hops, wall collisions, and disc collisions for two hard discs in a box, using a billiard model and ergodic theory, validated by simulations.

Contribution

It provides the first exact analytical expressions for collision and hop times in a two-disc system, connecting ergodic theory with molecular dynamics.

Findings

Analytical formulas match simulation results with high accuracy.

Exact mean times for hops, wall collisions, and disc collisions are obtained.

The approach links billiard models and ergodic theory to molecular dynamics.

Abstract

We study the molecular dynamics of two discs undergoing Newtonian ("inertial") dynamics, with elastic collisions in a rectangular box. Using a mapping to a billiard model and a key result from ergodic theory, we obtain exact, analytical expressions for the mean times between the following events: hops, i.e.~horizontal or vertical interchanges of the particles; wall collisions; and disc collisions. To do so, we calculate volumes and cross-sectional areas in the four-dimensional configuration space. We compare the analytical results against Monte Carlo and molecular dynamics simulations, with excellent agreement.

Equations

$$\begin{array}{rl}x:=\frac{x_1-x_2}{\sqrt{2}};& X:=\frac{x_1+x_2}{\sqrt{2}};\\ y:=\frac{y_1-y_2}{\sqrt{2}};& Y:=\frac{y_1+y_2}{\sqrt{2}}.\end{array}$$

$$⟨\tau ⟩=\frac{|Q|}{|A|}\cdot \frac{|S^{d-1}|}{|B^{d-1}|}\cdot \frac{1}{s}.$$

$$|S^{d-1}|=\frac{2{\pi }^{d/2}}{\Gamma (d/2)},$$

$$x_1=x_2\text{or}{y}_1=y_2.$$

$$V_{\text{free}}$$

$$=\int \mathrm{d}\mathbf{x}{1}_Z(\mathbf{x});$$

$$V_{\text{free}}=\underset{\begin{array}{c}x=-a\sqrt{2}\\ y=-b\sqrt{2}\end{array}}{\overset{\begin{array}{c}x=a\sqrt{2}\\ y=b\sqrt{2}\end{array}}{\iint }}\mathrm{d}x\mathrm{d}y\underset{\begin{array}{c}X=-a\sqrt{2}+|x|\\ Y=-b\sqrt{2}+|y|\end{array}}{\overset{\begin{array}{c}X=a\sqrt{2}-|x|\\ Y=b\sqrt{2}-|y|\end{array}}{\iint }}\mathrm{d}X\mathrm{d}Y{1}_{\{x^2+y^2\ge 2r^2\}}.$$

$$V_{\text{free}}$$

$$=16\underset{\begin{array}{c}x=0\\ y=0\end{array}}{\overset{\begin{array}{c}x=a\sqrt{2}\\ y=b\sqrt{2}\end{array}}{\iint }}\mathrm{d}x\mathrm{d}y\left(a\sqrt{2}-x\right)\left(b\sqrt{2}-y\right)1_{\{x^2+y^2\ge 2r^2\}},$$

$$I_1$$

$$=2a{b}^{2}r+\frac{1}{6}(a+b)(2r{)}^3-\frac{1}{32}(2r{)}^4-\frac{1}{4}\left(\pi ab+b^2\right)(2r{)}^2,$$

$$I_2$$

$$=\left(a^2-2ar+r^2\right)b^2.$$

$$V_{\text{free}}$$

$$=:V_{\text{prism}}-V_{\text{cyl}},$$

$$V_{\text{free}}=(w-2r{)}^2(h-2r{)}^2-\pi (w-2r)(h-2r)4r^2+\frac{32}{3}(w+h-2r)r^3-8^4.$$

$$\begin{gathered} V_{h/4<r<w/4}=32ab{r}^2\left[\arccos (b/2r)-\arccos (a/2r)\right] \\ +\frac{64r^3}{3}\left[a((b-a)/2r)-b(c/2r+\sqrt{4r^2-a^2}/2r)\right] \\ +16\left[a{b}^{2}c(4\sqrt{2}-1-\sqrt{2}/3)+c^2{b}^2(\sqrt{2}/3-1)\right] \\ -2r^2(b^2-a^2) \end{gathered}$$

$$A(S)=\int 1_Z(\mathbf{x})\delta (g(\mathbf{x}))\mathrm{d}\mathbf{x}.$$

$${\int }_{{\mathbf{R}}^d}f(\mathbf{x})\delta (g(\mathbf{x}))\mathrm{d}\mathbf{x}={\int }_{g^{-1}(0)}\frac{f(\mathbf{x})}{\Vert \nabla g(\mathbf{x})\Vert }\mathrm{d}S,$$

$$A_{\text{hop}}=\iiiint\limits_{\begin{subarray}{c}x_{1},x_{2}=-a\\ y_{1},y_{2}=-b\end{subarray}}^{\begin{subarray}{c}x_{1},x_{2}=a\\ y_{1},y_{2}=b\end{subarray}}\mathrm{d}{x}\,_{1}\mathrm{d}{x}\,_{2}\mathrm{d}{y}\,_{1}\mathrm{d}{y}\,_{2}\,\mathbf{1}_{\{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}\geq(2r)^{2}\}}\,\delta\left(\frac{y_{1}-y_{2}}{\sqrt{2}}\right).$$

$$A_{\text{hop}}=16\sqrt{2}b(a-r{)}^2.$$

$$g_{\text{coll}}(\mathbf{x}):=(x_1-x_2{)}^2+(y_1-y_2{)}^2-4r^2,$$

$$\nabla g_{\text{coll}}(\mathbf{x})=(2(x_1-x_2),-2(x_1-x_2),2(y_1-y_2),-2(y_1-y_2)),$$

$$A_{\text{coll}}$$

$$g_{\text{wall}}(\mathbf{x}):=x_1-a.$$

$$A_{\mathrm{wall}}=\underset{\begin{array}{c}{x}_2=-a\\ y_1,y_2=-b\end{array}}{\overset{\begin{array}{c}{x}_2=a\\ y_1,y_2=b\end{array}}{\iiint }}\mathrm{d}x{}_2\mathrm{d}y{}_1\mathrm{d}y{}_2{1}_{\{(a-x_2{)}^2+(y_1-y_2{)}^2\ge 4r^2\}}$$

$$A_{\mathrm{wall}}$$

$$A_{\text{walls}}$$

$$⟨{\tau }_{\text{hop}}⟩=\frac{3\pi }{4\sqrt{2}}\frac{2a^2{b}^2-2\pi ab{r}^2+\frac{a+b}{3}(2r{)}^3-r^4}{b\sqrt{2}(a-r{)}^2}.$$

$$⟨{\tau }_{\text{hop}}⟩\stackrel{r\to 0}{\sim }\frac{3\pi }{8\sqrt{2}}h.$$

$$⟨{\tau }_{\text{coll}}⟩=\frac{3\pi }{2\sqrt{2}}\frac{2a^2{b}^2-2\pi ab{r}^2+\frac{a+b}{3}(2r{)}^3-r^4}{2\pi abr-4(a+b)r^2+2r^3}.$$

$$⟨{\tau }_{\text{coll}}⟩\stackrel{r\to 0}{\sim }\frac{3}{8\sqrt{2}}\frac{wh}{r}.$$

$$⟨{\tau }_{\text{wall}}⟩=\frac{3\pi }{2\sqrt{2}}\frac{2a^2{b}^2-2\pi ab{r}^2+\frac{a+b}{3}(2r{)}^3-r^4}{a{b}^2-\pi /2b{r}^2+\frac{16}{3}{r}^3},$$

Full text

I Introduction

Collision times play a fundamental role in statistical mechanics, since properties such as mixing, reaction and diffusion rates depend on them Boltzmann (1872); Tolman (1938); VanKampen (2007). A paradigmatic model, introduced by Boltzmann almost 150 years ago Boltzmann (1872); Szasz (2000), consists of a fluid of hard spheres colliding with one other, either in a periodic torus or in a rectangular box.

The simplest such system that is non-trivial is that of two discs in a box. The discs move rectilinearly until they undergo elastic collisions, either with a wall of the box, or with one another; this has been called “inertial motion”, to distinguish it from Brownian motion Bowles et al. (2004). Phenomena such as transport, collision and reactions have been studied in such systems, both analytically Awazu (2001); Munakata and Hu (2002); Suh et al. (2005a) and numerically Suh et al. (2004, 2005b). The events of physical interest in this system are collisions of the particles, with the box or with one other, and hopping, in which two particles interchange position, either vertically or horizontally, and; this plays an important role in the dynamics of confined fluids, for example Bowles et al. (2004).

Previous work by Bowles et al. gives results for hopping for both inertial and Brownian dynamics Bowles et al. (2004), defined as the mean first-passage time, using arguments from transition state theory to provide expressions for the general behavior of these quantities, as a function of the disc radii. Further work has expanded along that line Suh et al. (2005a); Ball et al. (2009). There is also a statistical thermodynamic treatment by Munakata et al. , who study the partition function, pressure, and temperature in the system Munakata and Hu (2002); Uraganase and Munakata (2006).

In this paper we obtain exact analytical expressions, not for the first-passage time, but rather for mean inter-hop and inter-collision times of two discs under inertial (Newtonian) motion, as a function of the geometrical parameters of the system. Such a treatment is possible since the dynamics of $N$ hard discs in $d$ spatial dimensions is equivalent to a billiard model in the $Nd$-dimensional configuration space, i.e. a single point particle colliding elastically with suitable objects (“scatterers”) Szasz (2000). In particular, the system of two discs in a box is hyperbolic (chaotic) and ergodic when the discs are able to pass each other Simanyi (1999). Otherwise, the system decomposes into two or four disjoint ergodic components, on each of which the dynamics is still hyperbolic.

We can thus apply a key result from ergodic theory on mean collision times Chernov (1997) to obtain analytical expressions for times between hops, between wall collisions, and between disc collisions, in terms of volumes and areas in configuration space. We verify our results by comparing the geometrical calculations with rejection-sampling Monte Carlo simulations, and the dynamical results with molecular dynamics simulations. In certain asymptotic regimes, we recover previously obtained power-law exponents Bowles et al. (2004).

The main difficulty in the analysis is correctly accounting for several geometric factors with different origins that occur in both the analytical and numerical calculations. All numerical code, in the Julia language Bezanson et al. (2017), is available online rep .

II Model: Two hard discs in a rectangular box

We consider two hard discs with equal radius $r$, moving inside a box of width $w$ and height $h$; see Fig. 1. The discs move inertially in the absence of forces, following straight line trajectories, and undergo elastic collisions with one other and with the walls of the box.

We denote the position of the center of the $i$th disc by $\mathbf{x}_{i}=(x_{i},y_{i})$, and its velocity by $\mathbf{v}_{i}=(u_{i},v_{i})$ for $i=1,2$. Since the discs are hard, their centers are restricted to the region $(x_{i},y_{i})\in[-a,a]\times[-b,b]$, where $a:=a(r):=\frac{w}{2}-r$ and $b:=b(r):=\frac{h}{2}-r$.

The exclusion condition preventing the discs from overlapping is $(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}\geq(2r)^{2}$. It is thus useful to perform an orthogonal transformation (rotation) to the following new coordinates:

[TABLE]

In these coordinates, the configuration space is given by the following intervals: $x\in[-a\sqrt{2},+a\sqrt{2}]$ with $X\in[-a\sqrt{2}+|x|,a\sqrt{2}-|x|]$; and similarly for $y$ and $Y$, replacing $a$ by $b$. The non-overlapping constraint then becomes $x^{2}+y^{2}\geq 2r^{2}$. The horizontal and vertical coordinates transform independently from one other, and the Jacobian of this transformation is equal to $1$.

These constraints define a four-dimensional rectangular prism, in which is embedded an excluded cylinder with a three-dimensional surface (codimension 1). This cylinder has radius $r\sqrt{2}$ and lies on a diagonal between the $X$ and $Y$ axes. The prism surface is the outer boundary of the configuration space, while the cylinder is an excluded volume, the surface of which acts as a reflecting inner boundary. The dynamics of the two discs is equivalent to a billiard model in this 4-dimensional space, in which a point particle undergoes free flight until it hits a wall, where it undergoes an elastic reflection. The outer boundaries are flat, so the hyperbolicity is due to the inner semi-dispersing boundary Simanyi (1999), which corresponds to the collision of the two discs.

We take the mass of each disk as $m=1$, so that the kinetic energy is $\frac{1}{2}(\mathbf{v}_{1}^{2}+\mathbf{v}_{2}^{2})$. We restrict attention to the energy surface with $E=\frac{1}{2}$, so that the disc velocities satisfy $\mathbf{v}_{1}^{2}+\mathbf{v}_{2}^{2}=1$. Other values of the mass or energy correspond to a simple rescaling of the dynamics, with velocities differing by a factor of $\sqrt{2E/m}$, and corresponding factors in the times to be determined below.

II.1 Hopping dynamics

An example of hopping dynamics is shown in Fig. 2, for parameter values for which the two discs are trapped for long periods without being able to undergo horizontal hops. This was simulated using the algorithm described below. We plot the difference $x_{2}-x_{1}$ as a function of time $t$; horizontal hops correspond to the moments at which this quantity crosses through [math].

We see a bistable-type behavior corresponding to the two possible horizontal arrangements of the two discs. The goal of this paper is to calculate the mean time between adjacent hops. Note that, as seen in the figure, these statistics include pairs (and possibly longer sequences) of events that occur very close to each other in time, corresponding to consecutive hops without intermediate randomization. Although we do not necessarily wish to count these as true hopping events, it is not clear how to exclude them.

III Mean collision time for billiard models

A system of $N$ hard spheres confined by hard walls in a $d$-dimensional space may be treated as a billiard system in which a single point particle undergoes free motion between reflecting obstacles in a $(dN)$-dimensional configuration space Sinai (1970); Simanyi (1999); Chernov and Markarian (2006). This can be thought of as a mean return time to the $(d-1)$-dimensional (i.e. co-dimension $1$) cross-section given by the wall boundaries.

For a particle with unit velocity moving in a general billiard table, there is an exact expression for the mean return time between consecutive intersections of the particle with a given surface Chernov (1997):

[TABLE]

Here, $|Q|$ denotes the $d$-dimensional volume of the available space in the billiard; $|S^{d-1}|$ is the $(d-1)$-dimensional area of the unit sphere in $\mathbb{R}^{d}$, given by

[TABLE]

where $\Gamma(\cdot)$ is the gamma function; and $|B^{d-1}|$ is the volume of the unit ball in $\mathbb{R}^{d-1}$, given by $|B^{d-1}|=|S^{d-2}|/(d-1)$. $|A|$ denotes the $(d-1)$-dimensional cross-sectional area of the surface in question. The extra factor of $s$ is a “sidedness” factor, explained below. Note that Machta and Zwanzig Machta and Zwanzig (1983) used a similar method to derive an escape time across a virtual boundary by treating it as a recurrence time.

In our case, we are interested in the mean return time to several types of co-dimension-$1$ cross section. The first, giving the “hopping time”, is defined by the moment at which the two discs interchange their horizontal or vertical position, i.e. when they cross one of the two surfaces

[TABLE]

At such an event, the cross-section can be hit from either side, e.g. disk 1 can be traveling from left to right or the other way around at the moment when the interchange of positions take place in a horizontal hop. The area $A$ of the surface is, then, effectively twice as large, or equivalently the mean return time to the surface is reduced by half; this is taken into account by taking the “sidedness” factor $s=2$ in the formula for $\left\langle\tau\right\rangle$.

The other events of interest correspond to collisions of a specific disc with a specific wall, and the mutual collision between the two discs. In both of these cases, the collision surface can be reached only from one side in configuration space, corresponding to $s=1$.

To calculate the mean times of interest, it is thus necessary to calculate the 4-dimensional volume $V$ of the available space, and the 3-dimensional cross-sectional area $A$ for each event of interest.

IV Calculation of volume and areas

In this section, we calculate analytically the available volume and the cross-sectional areas required.

IV.1 Volume of available space

We denote by $Z:=\{\mathbf{x}\in\mathbb{R}^{4}:(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}\geq(2r)^{2}\}$ the complement of the cylinder in configuration space, where $\mathbf{x}:=(x_{1},x_{2},y_{1},y_{2})$. The four-dimensional available volume $V_{\text{free}}$ is then given by

[TABLE]

here, $\int\mathrm{d}\mathbf{x}$ denotes a four-dimensional integral over the whole volume, and $\mathbf{1}_{Z}$ is the indicator function of the set $Z$, given by $\mathbf{1}_{Z}(\mathbf{x})=1$ if $\mathbf{x}\in Z$, and $=0$ if $\mathbf{x}\notin Z$. Fig. 3 shows a diagram of the product of spaces that give rise to the whole configuration space. Recall that the dimensions $a$ and $b$ of the available configuration space are functions of the disc radius $r$, but we suppress this explicit dependence for simplicity.

It is easiest to represent the excluded cylinder in the coordinates defined in Eq. (1):

[TABLE]

Since $X$ and $Y$ do not appear in the integrand, the corresponding integrals are trivial, giving

[TABLE]

where we have used the symmetry visible in Fig. 3. Thus $V_{\text{free}}=16(I_{1}+I_{2})$, where $I_{1}$ is the region where the range of values of $y$ is affected by the exclusion condition, and $I_{2}$ is where the exclusion condition has no effect on that range. In figure 3 this would correspond to the dotted and white areas of the left diagram; see also Fig. 13 in the appendix. We have

[TABLE]

and

[TABLE]

Thus

[TABLE]

as previously obtained by Munakata and Hu Munakata and Hu (2002). It is useful to split this expression into the total volume of the prism, $V_{\text{prism}}=16a^{2}b^{2}$, and the cylindrical volume excluded by the overlapping condition, $V_{\text{cyl}}=16\pi abr^{2}-\textstyle\frac{64}{3}(a+b)r^{3}+8r^{4}$.

The substitutions $a(r)\leftarrow(w-2r)/2$ and $b(r)\leftarrow(h-2r)/2$ give us the volume as a function of the radius, for fixed table size:

[TABLE]

Note that the above formula is correct only when both vertical and horizontal hopping are possible, i.e., when $h,w>4r$. If this is not the case, the integration limits for $X$ and $Y$ in (7) are altered; see Appendix A for the corresponding results.

As an example, consider the case $w\geq h$, in which vertical hopping is possible but horizontal hopping is not. There are two sub-cases: if $h\leq w<2h$, there is a value of $r$ above which hopping is no longer possible, but the discs still fit in the table. For $2h\leq w$, vertical hopping is possible only up to $r=h/2$. The configuration space splits into disjoint components, but thanks to the symmetry of the problem, in some cases the cross-section areas and volumes become disjoint components sharing the same fraction of the total volume, making the transition continuous. In other cases, we instead end up with a discontinuity in the formulas, corresponding to a factor of 2 or 4.

For example, when $h/4<r<w/4$ defining $c=\sqrt{4r^{2}-b^{2}}$, we have

[TABLE]

When $r$ is larger than both $h/4$ and $w/4$, one must take into account a similar contribution which inverts the roles of $a$ and $b$.

We have verified our expressions for the available volume numerically using standard rejection-sampling Monte Carlo simulations: we generate uniform random positions for the disc centers in $[-a,a]\times[-b,b]$ and count the fraction of initial conditions for which the two discs do not overlap (rejecting those where is overlap). The results are shown in Fig. 4. We chose the values $w=1.5,h=1$ for the numerical integration since all the different cases of the volume formulas occur for these values: all hops possible, only vertical hops possible, no hops possible.

In the following, for simplicity we shall give results only for the case $r<\text{min}(\frac{w}{4},\frac{h}{4})$.

IV.2 Cross-sectional areas

The terms of the form $|A|$ in (2) are surface areas of 3-dimensional surfaces (manifolds) $S$ embedded in configuration space, defined by algebraic equations of the form $g(\mathbf{x})=0$, so that $S=g^{-1}(0)$ is the zero set of $g$. The surface area of $S$ is then given by

[TABLE]

To evaluate this, we use the following coarea formula (Hörmander, 1983, section 6.1)

[TABLE]

where the integral on the right-hand side is over the surface $g^{-1}(0)$. Zappa et al. (2018)

The appearance of the normalization factor $\|\mathbf{\nabla}g(\mathbf{x})\|$ can be understood intuitively by considering how to verify numerically these surface areas. One possible technique is to use rejection sampling to sample the volume given by $|g(\mathbf{x})|\leq\eta$ for some small value of $\eta$. Points will be accepted or rejected according to their orthogonal distance to the surface, which is the factor in the denominator in the above expression. Our case is simple in this respect: most of our surfaces are either parallel to the coordinates or are in some way “diagonal”, giving constant factors.

IV.2.1 Hopping

Vertical hopping occurs when the vertical positions of the two discs are equal: $y_{1}-y_{2}=0$. In the above formulation we define $g_{\mathrm{hop}}(x_{1},x_{2},y_{1},y_{2}):=y_{1}-y_{2}$, so that $\nabla g_{\mathrm{hop}}(\mathbf{x})=(0,0,1,-1)$ and $\|\nabla g_{\mathrm{hop}}(\mathbf{x})\|=\sqrt{2}$. Alternatively, we can integrate directly in the $(x,X,y,Y)$ coordinates, where $g(x,X,y,Y)=y$, and then return to the original coordinates. This results in the following expression:

[TABLE]

Carrying out the integrals (see Appendix A) gives

[TABLE]

As before, the formula is no longer valid for $r>w/4$. In this case, vertical hopping becomes impossible for a larger radius.

In Monte Carlo simulations, we count the proportion of successful placements of hard discs for which the distance $|y_{1}-y_{2}|$ is within a small tolerance of [math], but note that the factor $\sqrt{2}$ must also be taken into account in this calculation (see Appendix B). Results are shown in Fig. 5.

The results for horizontal hopping are obtained by interchanging $a$ and $b$.

IV.2.2 Disc collisions

The area which represents collisions between the two discs is the surface area of the cylinder that lies within the prism, given by

[TABLE]

so that

[TABLE]

with norm $\|\nabla g_{\text{coll}}(\mathbf{x})\|=2\sqrt{2}\sqrt{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}}$.

For $r<w/4$, we find for the corresponding area

[TABLE]

The numerics proceeds as before, checking which random configurations fall within a small tolerance from the collision condition, and plotting this as a fraction of the total volume. The result is shown in Fig. 6.

IV.2.3 Wall collisions

We restrict attention to collisions of disc $1$ with the right wall, for which

[TABLE]

The corresponding area, after taking the delta function into account, is

[TABLE]

The integration gives

[TABLE]

Once again, a simple Monte Carlo procedure verifies this result, shown in Fig. 7. This time, the correction factors of $\sqrt{2}$ do not appear, since the areas are orthogonal or parallel to the original position coordinates.

Taking into account the symmetry of the expression for either disc bouncing on either of the vertical walls, the area for this event is $4A_{\text{wall}}$. For horizontal wall collisions, the roles of $a$ and $b$ are switched. Thus the area for any impact on any wall is

[TABLE]

V Exact mean inter-event times

In this section, we use the results of the last section to give analytical expressions for mean inter-hop times, as well as mean wall collision and disc collision times. We tested the results against molecular dynamics simulations with uniform random non-overlapping initial conditions for the discs and standard billiard dynamics with elastic collisions, in which the trajectories of the disc centers are determined by their momentum and the first collision (either with a wall or between discs) is given by the event with the least positive collision time.

This dynamics produces a list of collision times with walls and between discs. Hopping times are inferred by locating consecutive collisions in which the sign of the difference of the relevant coordinates (horizontal or vertical) changes. The exact hopping time may then be determined using the momenta after the previous collision.

V.1 Mean hopping time

Inserting the results of the previous section into the formula for the mean times for crossing surfaces of section (2) gives exact mean inter-hop times. For vertical hops we have

[TABLE]

(Recall that here there is a factor $s=2$.)

In the limit of small disc radius, the discs have very infrequent interactions, and the result depends only on the table height:

[TABLE]

Also of interest is the limit when $r\to w/4$, at which point vertical hopping becomes impossible. The denominator goes to zero quadratically, while the available volume is still positive (except in the degenerate case $w=2h$). The exact expression is cumbersome, due to the fact that the general volume expression contains trigonometric functions and is unintuitive (see Section A.1), but the leading term in the numerator is $w^{2}(\frac{h}{2}-\frac{w}{4})^{2}$, and the denominator stays the same, so the asymptotic behavior is $\sim(r-\frac{w}{4})^{-2}$. The figure shows the comparison between numeric and analytic results for the complete valid range.

V.2 Mean disc collision time

For collisions between discs, we have, for the case $r<w/4$,

[TABLE]

As expected, this tends to infinity in the limit of small radius, with asymptotics

[TABLE]

For the case in which the discs narrowly fit inside the table we need to use the more cumbersome expression in (17) and the corresponding area. The time between collisions should go to zero; see figure 9. In this case, the function is smooth for the whole range of valid values for $r$.

V.3 Mean wall collision time

Lastly, for a specific disc colliding with a specific vertical wall we have

[TABLE]

in the case that $r<h/4$. The limiting form for small $r$ depends only on the width of the table and is exactly $3\pi w/2$. In the numerics this appears as the intersection with the $y$-axis.

When $r\approx h/4$ we have two limiting expressions since there is a discontinuity in the available volume; the difference is exactly a factor of two. In the comparison with numerics, we can see this jump in the function, accompanied by a larger error in the calculations (since it is more difficult to find correct starting positions). This is indicative of the system breaking up into smaller ergodic components.

VI Probability distributions

Since each type of event studied is a recurrence time to a certain surface in phase space, we expect to observe the standard exponential distribution for return times in chaotic systems Hirata et al. (1999). It is known that deviations from the exponential depend on the particularities of the system at hand Altmann and Kantz (2005).

As an example, we show the probability distribution of hopping times for different radii and find that the distribution depends both qualitatively and quantitatively on the radius, with the tails approximately exhibiting the expected exponential decay.

VII Conclusions

We have calculated exact analytic expressions for mean inter-event times for hopping and collision in the paradigmatic model of two discs in a hard box, by treating them as returns to suitable surfaces in configuration space, and applying a result from ergodic theory. The analytical results derived are in agreement with the limiting behavior obtained by Bowles et al.  Bowles et al. (2004) for first-passage events, and are in excellent agreement with molecular dynamics simulations.

The results may be extended to discs with different radii $r_{1}$ and $r_{2}$, leading to clumsier analytical expressions. It would also be possible to extend the results to calculate exact mean collision and hopping times for systems containing more discs, but the analytical calculations become more challenging; see e.g. Cao et al. (2004).

Appendix A General area and volume calculations

In this section, we find expressions for areas and volumes in all regimes. To simplify the formulas, we omit the dependence of $a$ and $b$ on $r$.

A.1 Volume

An implicit assumption was made on the limits of integration in (7). If $w,h>4r$, then the limits of integration are unaffected by the radius of the circles. In order to avoid the same positive term $16a^{2}b^{2}$, in the formulas, we work here with the excluded volume $V_{\text{cyl}}$, instead of the available one $V_{\text{free}}$.

We begin the derivation after integrating out $X$ and $Y$:

[TABLE]

A diagram helps us visualize the limits of integration. The most general case is (without loss of generality) $h<w<2h$; as the disc radius

• •

Horizontal and vertical hopping possible: $0<r\leq h/4$;

• •

Only vertical hopping possible: $h/4<r\leq w/4$;

• •

No hopping possible: $w/4<r<(h+w-\sqrt{2hw})/2$.

The largest possible radius is illustrated in Fig. 12.

We examine these regimes on the integration space. The first one is solved in the main text, but we repeat it here using a different coordinate system that proves useful. In Fig. 13 we present the three regimes as the shaded area where we perform the integration.

We evaluate the integral, carefully choosing the limits of integration. We perform the integral over the hatched area in polar coordinates $(\rho,\theta)$. We call the integral $V_{h}(r,\alpha,\beta)$ (where $h$ stands for “hatched”). As can be seen from figure 13, the regime where all hopping is possible corresponds to $\alpha=\pi/2,\beta=0$; the vertical hopping regime is where $\alpha<\pi/2,\beta=0$, and the no-hopping regime is where $0<\beta<\alpha<\pi/2$.

We rewrite the integral in terms of an angle $\theta$:

[TABLE]

The indicator function determines the limits for the $\rho$ integration variable. We set $\alpha,\beta$ as the other two limits and perform the integral over $\rho$, which does not change in the three regimes.

[TABLE]

Now we integrate over the $\theta$ variable:

[TABLE]

For the case in which all hopping is posible, the expression in (34) takes the values $\alpha=\pi/2,\beta=0$ and after multiplying by 16 both sides, we recover the Munakata and Hu formula. For the other two cases, we use that $\sin\alpha=b/r$ and $\cos\beta=a/r$.

Now we treat the dotted areas in the figures 13(b) and 13(c). Since they are triangular, they are more easily treated in Cartesian coordinates. We start with the upper triangular region. First, since we are inside the triangle the characteristic function becomes a simple integration limit:

[TABLE]

The lower dotted region in Fig. 13(c) gives a symmetric expression:

[TABLE]

These expressions account for all volume available in the configuration space, but cannot be used for calculation of all event times that interest us. We need also expressions that take into account that this space is divided into disjoint components, as some events become impossible in each of these subsystems. For example, (see the next section), disc 1 cannot hit the left wall if it started on the right and horizontal hopping is excluded. So we have to take into account that only half of the positions are available (due to symmetry) in (2).

Due to the symmetry of the problem, the available volume for each disjoint component of the dynamical system is equal; for example, if horizontal hopping is no longer possible ($w/4<r<h/4$), then there are two symmetric disjoint components: the system starting with disc 1 on the left and the one starting with disc 1 on the right. Both occupy the same phase space volume, as they are symmetric under interchange of labels. If $h<4\leq r$ then the system gets further divided into four disjoint components.

A.2 Area

The above procedure must be repeated for the calculation of areas for each surface of interest. We use a suitable Dirac delta to represent each codimension-1 collision event (contact of a disc with a wall or another disc). We multiply the characteristic function of the available space by the Dirac delta, and then again divide it into the three cases, namely, all hopping, only vertical hopping and no hopping possible, again referring to Fig. 13. Sometimes it turns out to be easier to obtain the integral over all configuration space and then exclude the part that corresponds to the overlapping condition. Again, the hatched area of the exclusion condition has a simpler representation in polar coordinates, and the triangular (dotted) regions can be treated in rectangular coordinates.

A.2.1 Hopping cross section

Fig. 14 shows the geometry for the hopping surface, given by $g_{\text{hop}}(\mathbf{x})=y_{1}-y_{2}=0$. We obtain

[TABLE]

The simplest approach to this integral is via the $y$ and $Y$ coordinates, since the surface is orthogonal to the $Y$ axis:

[TABLE]

Note that in (38) we do not use symmetry in $y$, due to the Dirac delta.

A.2.2 Disc collisions

For other cross-section areas we proceed in a similar manner. Disc collisions occur when $\sqrt{x^{2}+y^{2}}=r\,\sqrt{2}$, giving

[TABLE]

[TABLE]

We change to polar coordinates as we did in 32:

[TABLE]

giving

[TABLE]

We have used the same trick as in the previous subsection to obtain a general expression that works even when hopping is not possible. In the case of $r<w/4$ (all hopping possible) we have the substitution $\alpha=\pi/2$ and $b=0$, obtaining the result stated previously in eq. 22. It is practical to leave here the expression for 1/4 of the total area. The symmetry of the cases makes it so that when the phase space splits into disjoint components, the accessible volume and area scale in the same manner.

A.2.3 Wall collisions

For a collision with the wall, the $x_{i}$ give suitable coordinates. Once again, we have to be careful with the multiplicative factors that appear due to change of variables. Let us suppose that we want the area that represents contact of disc 2 with the right wall, so that $x_{2}-a=0$. Then

[TABLE]

We change variables as follows:

[TABLE]

Integrating over $x$, we obtain

[TABLE]

In the case that all hopping is possible there are no problems with the above derivations, but if vertical hopping is excluded then there are differences. If $r>h/4$, and disc 2 starts on the left, it will never be able to hit the right wall. So, actually, this is a subsystem of the whole system in which our general formula needs to be adjusted. When $h/4<r<w/4$, a disk cannot interchange left–right positions, but can still move across all vertical available space. After $r\geq w/4$, the system gets split again into two more disjoint components. So, when we calculate the available volume, we use only the part of the volume corresponding to the set where the event can occur.

Again, it is simpler to evaluate the last expression in eq. 53 in the excluded space and then to subtract that from the whole configuration space evaluation:

[TABLE]

The integrals are again routine to evaluate: the whole space part gives

[TABLE]

For the excluded part, we again evaluate the circular region in polar coordinates and the triangular regions, if they exist, in Cartesian coordinates. From the circular sector we get (see Fig. 13 and the following explanation for $\alpha,\beta$ variables):

[TABLE]

The upper triangle for $\alpha<\pi/2$ gives

[TABLE]

and the lower triangle a symmetric expression with $a$ instead of $b$.

Appendix B Numerical method for calculating surface area

All code used for the simulations is available, in keeping with the philosophy of open science.

For simplicity, we used traditional rejection sampling to calculate the free volume and cross-sectional areas.

The cross-sectional areas are calculated by fattening the surface $g(\mathbf{x})=0$ to $S_{\epsilon}:=\{\mathbf{x}:0\leq g(\mathbf{x})\leq\epsilon\}$, or in the case of the hopping cross-section $S_{2\epsilon}:=\{\mathbf{x}:\|g(\mathbf{x})\|\leq\epsilon\}$.

Thus

[TABLE]

is calculated by sampling points from the whole volume and rejecting those lying outside $S_{\epsilon}$.

As in the analytical calculation, the width of this strip must be calculated correctly by dividing by $\|\nabla g(x)\|$.

Acknowledgements

DPS thanks Sidney Redner for posing the question that led to this work, and the Marcos Moshinsky Foundation for financial support via a fellowship. Support is also acknowledged from grants CONACYT-Mexico CB-101246 and CB-101997, and DGAPA-UNAM PAPIIT IN117117.