Infinite horizon billiards: Transport at the border between Gauss and L\'evy universality classes

TL;DR

This paper analyzes transport in infinite horizon billiards, revealing non-Gaussian spreading behaviors, power-law tails, and the limitations of renewal assumptions, with implications for understanding Le9vy and Gauss universality classes.

Contribution

It provides analytical descriptions of particle spreading in billiards at the border between Le9vy and Gauss universality, highlighting the role of correlations and non-analytical collision time distributions.

Findings

Particles form power-law tails in density

Renewal assumption works for Lorentz gas but not for stadium channel

Slow convergence rate observed in particle spreading

Abstract

We consider transport in two billiard models, the infinite horizon Lorentz gas and the stadium channel, presenting analytical results for the spreading packet of particles. We first obtain the cumulative distribution function of traveling times between collisions, which exhibits non-analytical behavior. Using a renewal assumption and the L\'evy walk model, we obtain the particles' probability density. For the Lorentz gas, it shows a distinguished difference when compared with the known Gaussian propagator, as the latter is valid only for extremely long times. In particular, we show plumes of particles spreading along the infinite corridors, creating power-law tails of the density. We demonstrate the slow convergence rate via summation of independent identically distributed random variables on the border between L\'evy and Gauss laws. The renewal assumption works well for the Lorentz gas with intermediately sized scattering centers, but fails for the stadium channel due to strong temporal correlations. Our analytical results are supported with numerical samplings.

Equations

$$\underset{\tau \to \infty }{\lim }{\tau }^3\psi (\tau )={\tau }_0^2,$$

$$x=\sum _{n=1}^N{\chi }_n,$$

$$⟨\exp (i\bar{k}\bar{x})⟩=\exp \left\{N\ln \left[\tilde{f}\left(\frac{\bar{k}}{\sqrt{{\chi }_0^{2}N\Omega (N)/2}}\right)\right]\right\},$$

$$\tilde{f}(k)\simeq 1+\frac{1}{2}{\left({\chi }_{0}k\right)}^2\ln \left[{\left(C_f{\chi }_{0}k\right)}^2\right],$$

$$C_f=\exp \{\gamma$$

$$-{\int }_{{\chi }_0}^{\infty }\text{d}\chi \left[f(\chi ){\left(\frac{\chi }{{\chi }_0}\right)}^2-\frac{1}{\chi }\right]\},$$

$$⟨\exp (i\bar{k}\bar{x})⟩\simeq \exp \left\{\frac{{\bar{k}}^2}{\Omega (N)}\ln \left[\frac{2C_f^2{\bar{k}}^2}{N\Omega (N)}\right]\right\}.$$

$$\ln \left[\frac{N\Omega (N)}{2C_f^2}\right]=\Omega (N),$$

$$\Omega (N)=\left|{\text{W}}_{-1}\left(-\frac{2C_f^2}{N}\right)\right|.$$

$$\left|{\text{W}}_{-1}(\eta )\right|=L_1+L_2+\frac{L_2}{L_1}+\mathcal{O}\left(\frac{L_2^2}{L_1^2}\right),$$

$$⟨\exp (i\bar{k}\bar{x})⟩\simeq e^{-{\bar{k}}^2}\left[1+\frac{{\bar{k}}^2}{\Omega (N)}\ln \left({\bar{k}}^2\right)\right].$$

$$P(x,N)\simeq \frac{1}{\sqrt{\pi {\xi }^2(N)}}\exp \left[-\frac{x^2}{{\xi }^2(N)}\right]\times$$

$$\{1+\frac{1}{\Omega (N)}\{\left[2-\gamma -\ln (4)\right]\left[\frac{1}{2}-\frac{x^2}{{\xi }^2(N)}\right]$$

$$-\frac{1}{2}{\text{M}}^{(1,0,0)}\left[-1;\frac{1}{2};\frac{x^2}{{\xi }^2(N)}\right]\}\},$$

$$\mathbf{r}(t)-\mathbf{r}(0)=\sum _{n=1}^N{\mathbf{v}}_{n-1}{\tau }_n+{\mathbf{v}}_N{\tau }_b,$$

$${\Pi }_d\left(\mathbf{k},u\right)=\int {\text{d}}^{d}r{\int }_0^{\infty }\text{d}t{P}_d\left(\mathbf{r},t\right)e^{-ut+i\mathbf{k}\cdot \mathbf{x}}.$$

$${\Pi }_d\left(\mathbf{k},u\right)=⟨\frac{1-\hat{\psi }\left(u-i\mathbf{k}\cdot \mathbf{v}\right)}{u-i\mathbf{k}\cdot \mathbf{v}}⟩\frac{1}{1-⟨\hat{\psi }\left(u-i\mathbf{k}\cdot \mathbf{v}\right)⟩},$$

$$\hat{\psi }(u)={\int }_0^{\infty }\text{d}\tau \psi (\tau )e^{-u\tau },$$

$$⟨\cdots ⟩=\int {\text{d}}^{d}v\cdots F_d(\mathbf{v}),$$

$$\hat{\psi }(u)\simeq 1-⟨\tau ⟩u-\frac{1}{2}{\left({\tau }_{0}u\right)}^2\ln \left(C_{\psi }{\tau }_{0}u\right),$$

$$C_{\psi }=\exp \{\gamma$$

$$-{\int }_{{\tau }_0}^{\infty }\text{d}\tau \left[\psi (\tau ){\left(\frac{\tau }{{\tau }_0}\right)}^2-\frac{1}{\tau }\right]\}.$$

$$F_1\left(v_x\right)=\frac{1}{2}\left[\delta \left(v_x-V\right)+\delta \left(v_x+V\right)\right],$$

$${\tilde{P}}_1({\kappa }_x,t)\simeq \frac{2e^{-{\kappa }_x^2}}{{\xi }_1(t)}\left[1+\frac{1}{{\Omega }_1(t)}{\kappa }_x^2\ln \left({\kappa }_x^2\right)\right],$$

$${\Omega }_d(t)=\left|{\text{W}}_{-1}\left(-4d{C}_{\psi }^2\frac{⟨\tau ⟩}{t}\right)\right|,$$

$${\xi }_d(t)=\sqrt{{\tau }_0^2⟨v^2⟩\frac{t{\Omega }_d(t)}{d⟨\tau ⟩}},$$

$$P_d(\mathbf{r},t)=\int \frac{{\text{d}}^d\kappa }{(2\pi {)}^d}{\tilde{P}}_d\left(\mathbf{\kappa },t\right)\cos \left[\frac{2(\mathbf{\kappa }\cdot \mathbf{r})}{{\xi }_d(t)}\right].$$

$$P_1(x,t)$$

$$\times \{1+\frac{1}{{\Omega }_1(t)}\{\left[2-\gamma -\ln (4)\right]\left[\frac{1}{2}-\frac{x^2}{{\xi }_1^2(t)}\right]$$

$$-\frac{1}{2}{\text{M}}^{(1,0,0)}\left[-1;\frac{1}{2};\frac{x^2}{{\xi }_1^2(t)}\right]\}\}.$$

$$F_2\left(\mathbf{v}\right)$$

$$+\delta \left(v_x\right)\left[\delta \left(v_y-V\right)+\delta \left(v_y+V\right)\right]\}$$

$$+\frac{q}{4}\{\left[\delta \left(v_x-\frac{V}{\sqrt{2}}\right)+\delta \left(v_x+\frac{V}{\sqrt{2}}\right)\right]$$

$$\times \left[\delta \left(v_y-\frac{V}{\sqrt{2}}\right)+\delta \left(v_y+\frac{V}{\sqrt{2}}\right)\right]\},$$

$${\tilde{P}}_2(\mathbf{\kappa },t)$$

Full text

I Introduction

The infinite horizon Lorentz gas Lorentz is a paradigmatic model of deterministic classical transport, thoroughly studied by physicists Bouchaud ; Zaslavsky ; Cristadoro1 ; Cristadoro2 ; Pedro ; Klages ; Carlos and mathematicians Bunimovich1 ; Boldrighini ; Bunimovich2 ; Bleher ; Dahlqvist ; Szasz ; Zaharescu ; Chernov1 ; Chernov2 ; Chernov3 ; Marklof1 ; Dettmann1 ; Dettmann2 ; Marklof2 . It consists of an infinite periodic lattice of convex obstacles, and point-like particles which undergo elastic collisions with them. The most common configuration of the Lorentz gas model is composed of circular scatterers arranged into an infinite square lattice of unit spacing. Several important properties of this model, such as ergodicity Szasz and algebraic decay of the velocity correlations in time Bunimovich2 , have been rigorously proven. Importantly, Bleher Bleher showed that a particle’s position vector $\boldsymbol{r}(t)$ has a limiting Gaussian distribution when normalized correctly. More accurately, the quantity $\lim_{t\to\infty}[\boldsymbol{r}(t)-\boldsymbol{r}(0)]/\sqrt{t\ln(t)}$ is a two-dimensional Gaussian variable with zero mean and a covariance matrix which depends on the arrangement of scatterers. However, this asymptotic form is valid only when $\ln[\ln(N)]/\ln(N)=\epsilon\ll 1$, where $N$ is the number of collisions. To satisfy this condition, $N$ has to be extremely large (e.g., $\epsilon=0.01$ requires that $N\approx 10^{281}$). Bleher’s time is too large to be physically relevant Zaslavsky , while microscopic inter-collision times cannot describe transport processes happening on much larger time scales. As such, key features of the Lorentz gas model can only be seen for intermediate times that are most relevant for transport regimes. Unfortunately, currently there are no analytical results for these mesoscopic times.

In billiard models, an important characteristic of the transport is a presence (or an absence) of infinite horizons, corridors along which collision-less ballistic trajectories propagate, see Fig. 1. The packet of spreading particles in an infinite horizon billiard models exhibits two main features: the center part of the packet is approximately Gaussian, and the far tails are described by plumes of particles spreading along the infinite corridors. For a specific configuration of the Lorentz gas, we found in our previous Rapid Communication PRE2018 that the scatterers’ geometry is embedded in the cross-like shape of the spreading packet, see Figs. 1 (a) and 2. Here we wish to extend our theory to other configurations, showing its generality. For this aim, we consider the Lorentz gas with corridors forming a British flag-like structure for the packet of spreading particles, seen in Figs. 1 (b) and 3, together with a quasi one-dimensional transport in a chain of concatenated Bunimovich billiard stadiums Bunimovich1 ; Zaslavsky , seen in Figs. 1 (c) and 4. While the models are distinctive, along the corridors the far tails of the density decay spatially with a universal power-law, a feature well described by the Lévy walk model Levy1 ; Levy2 ; Levy3 ; Levy4 .

Our analysis is composed of two main ingredients, one of them is obtaining the aforementioned billiard systems’ full distribution of inter-collision times. For the Lorentz gas, results in the limits of $R\to 1$ and $\tau\to\infty$ (where $R$ is a scatterer’s radius and $\tau$ is an inter-collision time) were found by Bouchaud and Le Doussal Bouchaud . More recently, an asymptotic form in the limit of $R\to 0$ was established Dahlqvist ; Zaharescu ; Marklof1 . However, this important aspect in the characterization of transport is hardly discussed in the literature when finite sized scatterers are considered. This distribution’s probability density function (PDF) exhibits a far tail obeying

[TABLE]

which is valid for both limits of large and small scatterers. Therefore, we use a Lévy walk model with exponent $-3$ for the other ingredient, which is calculating the density of spreading particles. A principal issue here is that in a Lévy walk approximation one uses a renewal assumption, i.e. one neglects correlations between consecutive collisions. It was shown that for the cross-like configuration of the Lorentz gas and in the limit of large scattering centers, this condition is nullified as an effective trapping mechanism emerges Cristadoro1 . In contrast, when the scattering centers are not too large, the Lévy walk with the obtained cumulative distribution function (CDF) of the waiting times works perfectly, as demonstrated below. Deviations from the renewal theory do exist for the stadium channel model where correlations are strong, as we show below. However, this does not imply that the Lévy walk model is not predictive here as well, and in fact we find the opposite. Another primary issue is caused since these two systems are operating on the border between Gauss and Lévy central limit theorems, due to the exponent $-3$ which causes a logarithmic divergence of $\psi(\tau)$’s second moment. Thus, an ultra-slow convergence rate problem arises that can be understood via a toy model: summation of independent identically distributed (IID) random variables (RVs), with a common symmetric PDF that decays algebraically with a power of $-3$ for large argument. Already here we encounter the same type of convergence problem discussed in the first paragraph, which is neutralized using what we call Lambert scaling approach. This is a crucial step for these systems, as it allows us to compare finite time simulations with our theory for duration regimes where previous results do not hold.

The rest of this paper is organized as follows: in section II we provide an example of the Lambert scaling for sums of IID RVs. In section III we derive the main formula of the spatial PDFs by using the Lévy walk model. In section IV we obtain the CDFs of inter-collision times for the Lorentz gas and stadium channel. We discuss our results in section V.

II A simplified case of Lambert scaling

We now consider the problem of summation of IID RVs, drawn from a power-law distribution. We work at the border between Gauss and Lévy central limit theorems, which is clearly related to the exponent $-3$ in Eq. (1). Some aspects of this by far simpler approach are important for the discussed billiard models. In particular, at this transition we find a critical slowing down in the sense that convergence to the Gaussian limit theorem is ultra-slow Greengard , a problem which is resolved by Lambert scaling. Consider a sum of $N\gg 1$ IID RVs

[TABLE]

where the summands are drawn from a common symmetric PDF which obeys $f(\chi\to\infty)\simeq\chi_{0}^{2}/\chi^{3}$. We define the scaled sum as $\bar{x}=x/\sqrt{\chi_{0}^{2}N\Omega(N)/2}$ with $N\Omega(N)$ being a scaling parameter, soon to be determined. We use the characteristic function

[TABLE]

where $\tilde{f}(k)$ is the Fourier transform of $f(\chi)$. Assuming that $\Omega(N)$ monotonically increases with $N$, the small $k$ behavior of $\tilde{f}(k)$ is considered

[TABLE]

which is derived in appendix A. The first term is the normalization, while the second is related to the power-law tail of $f(\chi)$, with $C_{f}$ being

[TABLE]

where $\gamma\approx 0.5772$ is Euler’s constant. Inserting Eq. (4) into Eq. (3) and expanding, we get

[TABLE]

We now determine the slowly increasing function $\Omega(N)$ with the choice

[TABLE]

which yields

[TABLE]

Here $\text{W}_{-1}(\eta)$ is the secondary branch of the Lambert W-function Functions , defined for $\eta\in[-1/e,0)$ by the identity $\text{W}_{-1}(\eta)=\ln[\eta/\text{W}_{-1}(\eta)]$. The Lambert function has the following expansion as $\eta\rightarrow 0^{-}$

[TABLE]

where $L_{1}=\ln(1/|\eta|)$ and $L_{2}=\ln[\ln(1/|\eta|)]$. Equations (8) and (9) yield $N\Omega(N)\simeq N\ln(N)$ when $N\to\infty$, reproducing the well-known Gnedenko-Kolmogorov scaling $\sqrt{N\ln(N)}$ Limit . However, for this to be of relevance, one must demand $\ln(N)\gg\ln[\ln(N)]$, which makes the convergence to this mathematical limit ultra-slow. Our Lambert scaling approach unites all of the $N$-dependent logarithmic terms into a single function, thus resolving this problem. We therefore expand Eq. (6) for large $\Omega(N)$, keeping terms up to sub-leading order during the calculation

[TABLE]

Inverting Eq. (10) back to position space, we reach our first result, the PDF of $x$

[TABLE]

with $\xi(N)=\sqrt{2\chi_{0}^{2}N\Omega(N)}$. Here $\text{M}(\cdots)$ is Kummer’s confluent hypergeometric function Functions , and the superscript over M denotes its derivative with respect to its first argument. Figure 5 shows a good match between Eq. (II) and the sum Eq. (2)’s PDF for two different common density functions, $f(\chi)$ and $g(\chi)$ which are defined in the caption, with $N=10^{4}$. The Kummer function behave asymptotically as $\text{M}^{(1,0,0)}(-1,1/2,\eta^{2})\simeq-\sqrt{\pi}\exp(\eta^{2})/\eta^{3}$, suggesting that $P(x,N)\simeq N\chi_{0}^{2}/x^{3}$. Thus, the PDF of the sum $x$ reproduces the same power-law tails as of the original common PDF $f(\chi)$. As the distribution of traveling times between collisions Eq. (1) exhibit the same heavy-tail exponent of $-3$, our next step is implementing Lambert scaling to the Lévy walk model.

III Lambert scaling of the Lévy walk model

The $d$-dimensional Lévy walk model Levy2 ; Levy4 is defined as follows: a random walker is placed at $\boldsymbol{r}(0)$ on time $t=0$. Its movement consists of segments of ballistic motion with constant velocity, separated by collision-like events which induce a change in the velocity’s magnitude and/or direction. The process lasts for a fixed duration, which is the measurement time $t$. The model employs two PDFs in order to determine the displacement during each of the ballistic motion epochs. The velocity of each segment is drawn from a PDF $F_{d}(\boldsymbol{v})$, whose moments are all finite, and is further assumed to be symmetric with respect to each of the components $v_{j}$ where $1\leq j\leq d$ (such that its odd moments vanish). The time duration of each ballistic section is drawn from a PDF $\psi(\tau)$. The movement continues until the allotted measurement time is met, thus the number of collisions $N$ in $[0,t]$ is random. This yields the total displacement as

[TABLE]

where $\tau_{n}$ is the traveling time of the $n$th walking epoch, $\boldsymbol{v}_{n}$ is the velocity after the $n$th collision, and the initial conditions $\boldsymbol{r}(0)$ and $\boldsymbol{v}_{0}$ are randomly chosen. The traveling times and velocities $\{\tau_{n},\boldsymbol{v}_{n}\}$ (with $1\leq n\leq N$) are IID RVs, with the last movement duration being $\tau_{b}=t-\sum_{n=1}^{N}\tau_{n}$. Notice that the measurement time divided by the mean time between collisions, $t/\langle\tau\rangle$, and the lengths of inter-collision travel, $\{\boldsymbol{v}_{n-1}\tau_{n}\}$, roughly correspond to $N$ and $\{\chi_{n}\}$ from the previous section, respectively. We denote the particle’s speed with $V$, which is kept unchanged in the billiards systems due to the collisions’ elasticity (in the numerical simulations $V=1$). Let us denote the probability to find the walker at position $\boldsymbol{r}$ on time $t$ as $P_{d}(\boldsymbol{r},t)$, and let $\Pi_{d}(\boldsymbol{k},u)$ be its Fourier and Laplace transform

[TABLE]

An exact expression of $\Pi_{d}(\boldsymbol{k},u)$ is given by the Montroll-Weiss equation MW

[TABLE]

where $\hat{\psi}(u)$ is the Laplace transform of $\psi(\tau)$, defined as

[TABLE]

and

[TABLE]

with the above integral carried over all velocity space.

We now direct the reader’s attention to two points before advancing with the calculation. First, notice that the assumption of finite moments for the velocity distribution $F_{d}(\boldsymbol{v})$ is crucial to our theory. For example, a diverging second moment yields Lévy statistics in the bulk instead of a Gaussian, whereas for the Lorentz gas it is rigorously proven to be the latter case Bleher . Second, in Ref. Levy4 the authors discuss two different variations of Lévy walk: the velocity model and the jump model. For the former, particles move constantly until the measurement time ends, their last traveling epoch being $\tau_{b}$. However, for the latter particles are missing this final segment of walk. Therefore, the two models differ only when considering approximation theories of the far tail (under the condition of finite mean time between collisions). As our theory is intended to approximate the bulk of the PDF $P(\boldsymbol{r},t)$, these two aforementioned cases are indistinguishable. Technically speaking, this difference is manifested by a different numerator of Eq. (14), which does not change its approximated versions which appear below, Eqs. (20) and (III). Further, in the Lorentz gas a particle’s velocity is unity at any moment of travel, therefore the velocity model is more suitable if compared to the jump model.

To approximate Eq. (14) we use the asymptotic behavior of $\psi(\tau)$ Eq. (1), which implies that for small $u$

[TABLE]

which is derived in appendix B. Here, the first term is the normalization, $\langle\tau\rangle$ is the mean time between collisions, and the last term is related to the power-law tail of $\psi(\tau)$, with $C_{\psi}$ being

[TABLE]

We assume a scaling of $u\sim k^{2}L(k)$, where $L(\cdots)$ is some logarithmic-like function and $k=|\boldsymbol{k}|$. This suggests that $u\ll k$ when $k\rightarrow 0$, and as such we can expand the Montroll-Weiss Eq. (14) in the small parameter $u/kv$ where $v=|\boldsymbol{v}|$, see appendix C. For the one-dimensional stadium channel model, we apply the following distribution of velocities

[TABLE]

where in our case $V=1$. We obtain in appendix C the following Fourier expansion for the distribution

[TABLE]

with

[TABLE]

where here and for the following Lorentz gas $\langle v^{2}\rangle=V^{2}$. The walker’s position PDF in $d$ spatial dimensions is given by the inverse Fourier transform

[TABLE]

Inserting Eq. (20) into Eq. (22) results with

[TABLE]

The subscript $1$ stands for one dimension, namely this result should hold for the stadium channel when we coarse grain over the channel’s width. We see some similarities with the problem of summation of IID RVs. For example, Kummer’s function appears in both problems as a correction to the leading term. There is however a major difference between the two cases: here $P(\boldsymbol{r},t)=0$ for $\boldsymbol{r}>t$, which is different from the problem of summation of IID RVs. This is clearly due to the particles’ finite speed. For the two-dimensional Lorentz gas model with two/four infinite horizons we have

[TABLE]

where $q\geq 0$ is a parameter to be determined later, which encodes the probability of a particle to be found in the far tail of the diagonal corridors (for two open horizons, a cross-shape, one has $q=0$). We get from Eq. (14)

[TABLE]

We use a $\pi/4$ rotation transformation $\kappa_{\pm}=(\kappa_{x}\pm\kappa_{y})/\sqrt{2}$ to calculate the integrals over Eq. (III)’s last two terms, and obtain

[TABLE]

where the subscript $2$ stands for two dimensions. This solution is presented in Figs. 2 and 3, with the relevant parameters, namely $\tau_{0}$, $\langle\tau\rangle$, $C_{\psi}$ and $q$, obtained from $\psi(\tau)$ [more precisely we extract them out of $\text{CDF}(\tau)$]. Thus, we continue with deriving exact expressions for the distribution of traveling times $\psi(\tau)$ for the Lorentz gas model with two/four open infinite corridors and for the stadium channel model. To refrain from cumbersome formulas, we omit some of the next section’s derivations. Please refer to appendix D for more details.

IV The distribution of traveling times

IV.1 Lorentz gas model

To obtain the distribution of traveling times we define a $2d$ cubic lattice of constant $1$ occupied with circular scatterers of radius $R$, such that the center of each circle is located at a grid point, see Fig. 1 (a) and (b). Let the lattice coordinates of each scatterer be $(n,m)$, where $n$ and $m$ are integers. We focus on the origin, and assume that the particle has just collided with the $(0,0)$ scatterer. We define the collision’s impact parameter and recoil direction as $b$ and $\beta$ respectively, see Fig. 6 (a), where the ranges of values for these two parameters are $[-R,R]$ and $[0,2\pi)$, respectively. We denote as $\tau_{n,m}^{*}(b,\beta,R)$ the time duration until the following collision, and since $V=1$ it is also the distance traveled. Here, we assume that the next scatterer to be collided with is the $(n,m)$ one. The expression obtained for $\tau_{n,m}^{*}(b,\beta,R)$ is

[TABLE]

We then write the PDF $\psi(\tau)$ as

[TABLE]

where the factors of $1/2\pi\times 1/2R$ are the distributions of $\beta$ and $b$ respectively, which are both uniform due to ergodicity Bouchaud . In Eq. (28), the summation is carried over all relevant integers (see Fig. 7), and the integration boundaries (IBs) need to be found. Therefore, we define the discriminant of Eq. (IV.1)

[TABLE]

Our starting point for obtaining the IBs is to notice that the rooted quantities in Eq. (IV.1) must be positive. Actually, $R^{2}-b^{2}\geq 0$ means that the particle path indeed intersects the $(0,0)$ circle, which is an initial assumption here, and a positive discriminant in Eq. (29) means that the particle does collide with the $(n,m)$ scatterer. However, in order for the particle to reach the $(n,m)$ scatterer, it must as well not collide with another circle along its path. We therefore use a system of inequalities which are drawn from Eq. (29) to determine $\beta$ and $b$’s IBs. Assuming that $1/\sqrt{8}<R<1/2$, one has a single pair of infinite corridors, and the only scatterers which are reachable to the particle are those with lattice indexes that obey $n=1$ or $m=1$, see the textured colored circles in Fig. 7 (a). Symmetry considerations allow us to break the problem into eight areas composed of three components each, and we choose to focus on $m=1$ and $n\geq 0$, see the black rectangle-encircled area in Fig. 7 (a). This area’s first component is the nearest neighbor scatterer $(0,1)$. Here, Eq. (29) suggests that $b$ obeys

[TABLE]

The above condition ensures that the particle does collide with the scatterer $(0,1)$ on the next collision. In addition, this domain should be intersected with $|b|\leq R$, which ensures us that the particle indeed originated from the $(0,0)$ scatterer. Similarly, for the second component, the next to nearest neighbor scatterer $(1,1)$, the system of inequalities which stems from Eq. (29) reads

[TABLE]

and of course $|b|\leq R$ as before. This ensures that the particle collides with the $(1,1)$ scatterer (first condition), but not with the $(1,0)$ or the $(0,1)$ circles (second and third conditions, respectively) that will otherwise block its path. Lastly, the third component in this area is the distant neighbors cluster $n>1$, which is dominated by the inequalities

[TABLE]

as well as $|b|\leq R$. This ensures that the particle collides with the $(n,1)$ scatterer, but avoids from the $(1,0)$ or the $(n-1,1)$ ones, which can block its path. Analyzing these inequalities in appendix D, we obtain

[TABLE]

for the lower $\beta$ IB,

[TABLE]

for the upper $\beta$ IB, and

[TABLE]

for the IBs of $b$, where

[TABLE]

Equations (33-36) together with Eq. (28) provide a full description of $\psi(\tau)$, and alternatively by an integration over $\tau$, a full description of $\text{CDF}(\tau)$. Figure 8 (a) depicts the CDF obtained from these equations for $R=0.4$ with its respective numerical simulations counterpart, where an excellent match can be witnessed. With these equations in mind, we derive an exact expression for $\tau_{0}$, and numerical values of $\langle\tau\rangle$ and $C_{\psi}$ for $R=0.4$, see appendix E. For $\tau_{0}$ we obtain

[TABLE]

This matches the limiting result obtained in Ref. Bouchaud for $R\to 1/2$. It is also in perfect match to our previous exact result for $\tau_{0}$ which we obtained in Ref. PRE2018 using a different indirect approach. For the aforementioned numerical values, we get $\langle\tau\rangle\approx 0.62155$ and $C_{\psi}\approx 4.4802\times 10^{-4}$, where the former has a relative error of $0.021\%$ to the known rigorous result

[TABLE]

which is mentioned in Dettmann1 . These values provide excellent results for the numerical simulations of the position’s PDF when used as an input for the Lévy walk approximation, as seen in Figs. 2 and 9.

Next we compute the case of four open horizons for the Lorentz gas model, see Figs. 3 and 10. Assuming that $\sqrt{20}\leq R<1/\sqrt{8}$, there are two additional open horizons, namely the two main diagonals. It turns out that the formulas we obtained above for the $m=1$ and $n\geq 0$ stripe are valid here as well, excluding the $(2,1)$ scatterer. Since $R$ is now smaller, possible trajectories emerge for the diagonal directions, and we choose to focus on movements which end at the $(m+1,m)$ circles where $m\geq 1$, see the ellipse-encircled area in Fig. 7 (b). Thus, we need to adjust the $(2,1)$ scatterer’s upper IB, as it is now counted in two distinct sets of scatterers. The specifics are detailed in appendix D, and we obtain the correction

[TABLE]

The diagonal area’s inequalities are obtained yet again from Eq. (29). For the diagonal part of $(2,1)$, we have

[TABLE]

This circle is similar to the $(0,1)$ one in the previous case, as there are no possible obstacles for the (diagonal) nearest neighbor. The rest of the IBs are found by analyzing

[TABLE]

in a similar way as was done previously, where $m>1$. This time, the possible scatterers to block the particle’s path are $(m,m-1)$ and $(1,1)$. We obtain the following

[TABLE]

for the upper $\beta$ IB,

[TABLE]

for the lower $\beta$ IB, and

[TABLE]

for the IBs of $b$, where

[TABLE]

Using Eqs. (IV.1-45), we calculate the $\text{CDF}(\tau)$ for $R=0.3$, which is plotted in Fig. 8 (b), where an excellent match to the simulations can be seen. We derive from these equations an exact expression for $\tau_{0}$, obtaining

[TABLE]

Equation (46) is used to define $q$, which determines the relevant weight of the velocities’ PDF in the Lévy walk formalism, along respective directions of the infinite corridors. Recalling that $q$ determines the probability of a particle to be at the diagonal corridors, we define this parameter as the ratio between the diagonal corridors’ contribution for the behavior $\psi(\tau)\sim\tau^{-3}$ and the overall $\tau_{0}^{2}$, and find that

[TABLE]

when $1/\sqrt{20}\leq R<1/\sqrt{8}$. We also get from $\psi(\tau)$ numerical values for $\langle\tau\rangle$ and $C_{\psi}$ for this specific value of $R$, see appendix E. We find that $\langle\tau\rangle\approx 1.1947$, with a relative error of $0.059\%$ to the rigorous result Eq. (38), and also $C_{\psi}\approx 1.5250\times 10^{-2}$. These values provide excellent results for the numerical simulations of the position’s PDF when used as an input for the Lévy walk approximation, as seen in Figs. 3 and 10.

IV.2 Stadium channel model

Let us now precisely define the notation used for the stadium channel. As we have a lower and upper boundaries for this pipe structure, we define two parallel one-dimensional straight lattices of constant $2$ which are separated by a distance $D$. These are occupied with circular stadiums of radius $1$, such that the center of each stadium is located at a grid point, see Fig. 1 (c). Let $(2n,m)$ denote the center of a given stadium, where $n$ is an integer and $m$ can take two possible values, [math] and $D$ [see Fig. 6 (b)]. We focus on the origin, and assume that the particle has just been scattered from the lower wall’s $n=0$ stadium. We define the collision’s impact parameter and recoil direction as $a$ and $\alpha$ respectively, see Fig. 6 (b), where the ranges of values for these two parameters are $[-1,1]$ and $[0,2\pi)$, respectively. We denote as $\tau_{2n,m}^{*}(a,\alpha)$ the time duration until the following collision, and since $V=1$ it is also the distance traveled. Here, the next stadium to be collided with is $(2n,m)$, where the integer $n$ can be regarded as a semicircle’s numbering, while $m$ denotes the top or bottom wall. Notice that the particle is able to reach the stadium of origin, namely $n=m=0$. We obtained for $\tau_{2n,m}^{*}(a,\alpha)$

[TABLE]

The only semicircles which are reachable to the particle are those located at the upper row, or the origin stadium. Symmetry considerations allow us to break the problem into two areas, and we choose to focus on $n\geq 0$, see Fig. 1 (c). We then write the PDF $\psi(\tau)$ as

[TABLE]

where the factors of $1/2\pi\times 1/2$ are the distributions of $\alpha$ and $a$ respectively. As the corresponding parameters of the Lorentz gas $\beta$ and $b$ are known to be uniform Bouchaud , here we assume the same for $\alpha$ and $a$. The particle can travel into the $n\leq 0$ area, as well as from top to bottom, hence the multiplicative factor of $4$. To obtain the IBs, we use a similar scheme as for the Lorentz gas. However, there is a major difference between the two cases. Due to convexity, there are two possible points of origin/collision for the particle’s trajectory. In the Lorentz gas case, there was no need to differentiate between these two options during the calculations, thus we used the discriminant to tell whether a given scatterer was hit/missed. In the stadiums channel case, the discriminant is of no use as the stadiums are semicircles, and the discriminant cannot differentiate between a true stadium and a continuation of its wall to a complete circle. Therefore, here we define the vertical axis coordinate of the origin and target points, $y_{0}(a,\alpha)$ and $w_{2n,m}(a,\alpha)$ respectively

[TABLE]

and extract the needed inequalities from them instead. This time we demand that $w_{2n,m}(a,\alpha)-m$ will be non-negative for the upper stadiums, and non-positive for the origin stadium. This replaces the demand of a positive discriminant in the Lorentz gas case. We also demand that $y_{0}(a,\alpha)$ will be non-positive, which is analogous to $|b|\leq R$ in the Lorentz gas case. The first component of the chosen area is the origin semicircle, for which Eq. (IV.2) dictates that $a$ obeys

[TABLE]

The upper IB of $\alpha$ is set to $\pi/2$, and the lower to $-\pi/2$, which is possible due to symmetry. The second component is the $n=0$ upper stadium, for which we have

[TABLE]

The upper IB of $\alpha$ is set again to $\pi/2$ for symmetry considerations, and the lower is found by analyzing Eq. (52). Lastly, in order for the particle to reach the $n>0$ upper semicircles, it must not collide with another stadium wall along its path. For this third component, the conditions are

[TABLE]

The upper row in Eq. (IV.2) ensures that the particle originated from and arrived to the correct points, while the lower row prevents the top $2n-2$ and bottom $n=0$ semicircles from blocking the particle’s path. Analyzing these inequalities in appendix D, we obtain for the origin stadium as the target, $n=m=0$

[TABLE]

and for the top $n\geq 0$ and $m=D$ row of stadiums as targets

[TABLE]

for the lower $\alpha$ IB,

[TABLE]

for the upper $\alpha$ IB, and

[TABLE]

for the IBs of $a$, where

[TABLE]

Figure 11 depicts the CDF obtained from Eqs. (55-58) for $D=1$ with its respective numerical simulations counterpart, where an excellent match can be witnessed. We also use these to derive an exact expression for $\tau_{0}$, and the numerical values of $\langle\tau\rangle$ and $C_{\psi}$ for $D=1$, see appendix E. For $\tau_{0}$, we obtain

[TABLE]

We also find $\langle\tau\rangle\approx 2.57016$, with a relative error of $0.005\%$ to the simulation result $\langle\tau\rangle\approx 2.57031$, and $C_{\psi}\approx 1.0903\times 10^{-5}$. However, these values do not provide a correct description for the Lévy walk approximation, due to the renewal assumption being nullified by strong temporal correlations discussed below. Nonetheless, one can fit the Lévy walk approximation to the simulations data using a two parameters fit ($C_{\psi}^{2}\langle\tau\rangle$ and $\langle\tau\rangle/\tau_{0}^{2}$), thus find “effective" constants. These turn out to describe the problem well, as seen in Fig. 4 (a). We have verified that the values obtained for the constants by fitting do not change over time, see Fig. 4 (b).

We would like to direct the reader’s attention to the different geometry of the Lorentz gas CDF and the stadium channel CDF. Both CDFs have qualitatively different shapes, being convex/concave for the Lorentz gas/stadium channel, see Figs. 8/11 respectively. This might be related to the scatterers’ shape in the two models, which is convex/concave for the Lorentz gas/stadium channel, respectively. We believe this phenomenon is general as the geometry of the scattering centers is clearly embedded in this basic distribution, however we leave this intriguing point for a future work.

V Discussion and summary

Returning to section II, we address an issue we previously disregarded with the Lambert scaling approach. The reader may have noticed that the choice of the scaling function $\Omega(t)$ is not unique, but can be determined up to a constant. More accurately, one may choose to separate the logarithmic term into two at an arbitrary point, as one can always write

[TABLE]

Recall that we used the first term of Eq. (60)’s right hand side with $\eta=1$ to derive the Lambert scaling, see Eq. (7), while expanding the $\bar{k}^{2}$-containing exponential term, see Eq. (10). Here $\eta$ is a free parameter which cannot be determined uniquely by the aforementioned steps alone. We found that taking $\eta=1$ produces good results for $P(x,N)$ and alternatively for $P_{d}(\boldsymbol{r},t)$, see Figs 4, 5, 9, and 10. Of course, if one sums the complete asymptotic series Eq. (II) $\eta$ vanishes, but then the result diverges.

Looking back at Figs. 2 and 3, the reader may notice that the shape at $|\boldsymbol{r}|=r\sim t$ of the analytical results mismatches the simulations’. This can be explained via our scaling assumption, $u\sim k^{2}\ln(k)$. It suggests that $k\gg u$, namely our approximation is for displacement that obey $r\ll t$. Nonetheless, it holds well at the distribution’s infinite corridors, $y=0$ and $x=0$ for two open horizons and additionally $y=x$ and $y=-x$ for four. Utilizing the recently gaining attention infinite covariant density ICD could probably supply one with tools to approximate edge phenomena such as this.

Finally, we address the reader to interesting correlation patterns of the stadium channel model, seen in Fig. 12 (b). We plot points of consecutive traveling times $(\tau_{n+1},\tau_{n})$, for the Lorentz gas and the stadium channel, obtained from numerical samplings and analytical results. Both models exhibit a phenomenon of clear pointless areas on the graph, which is caused by the plateaus in the CDFs Figs. 8 and 11, which in turn correspond to a vanishing PDFs $\psi(\tau)$. Indeed, due to the discrete nature of the scattering centers, there exist certain durations of travel that are not possible (in a disordered system this non-analytical behavior would vanish). Technically speaking, large traveling times become semi-discrete, which is due to small parameter spaces $\{b,\beta\}$ [see Eqs. (33-35)] and $\{a,\alpha\}$ [see Eqs. (55-57)] which support these long trajectories. Interesting patterns emerge from roughly $\tau_{n}\leq 4.25$ for the stadium channel simulations, see Fig. 12 (b). These patterns suggest a highly correlative system. Indeed, when no correlations are present, one expects to find square-shaped patterns that correspond to independence of the axes, e.g. as seen for the Lorenz gas case (both analytical and numerical). However, the stadium channel’s simulations reveal a complete opposite. Take for example the line $\tau_{n+1}=\tau_{n}$, extending roughly up to $\tau_{n}\approx 1.5$. Examining the raw data, we found repeated instances of the form $\tau_{1},\tau_{1},\tau_{1},...$, where the number of elements depends on the size of $\tau_{1}$. Each instance always summed to a value of $\leq\pi/2$, the length of a semi-circle arc, see Fig. 13 (right cell). In these cases, a particle arrives at an almost vertical angle to one end of an arc and then propagates along it, similarly to gallery whispering modes of wave propagation. Other patterns correspond to different repeated instances, e.g. a particle performing “periodic orbits" [see Fig. 13 (left cell)] will have $\tau_{n}\approx 2+D=3$, which fits the middle rectangle. We computed the correlation functions for both models, and found that the Lorentz gas one decay much faster than the stadium channel one (not shown). Indeed, bounds on the temporal decay of the CDFs were obtained, showing that the Lorentz gas one decays with time in a stretched-exponential form Chernov1 , while the stadium channel one decays in a polynomial form Chernov2 . Thus, it is of no surprise that the renewal condition does not hold for the stadium channel.

In conclusion, we presented a geometrical method which yields the CDF of traveling times between collision events for billiard systems. We implemented this method for the Lorentz gas with infinite horizon and the stadium channel models. The approach is based on a shadowing effect (as some scatterers cannot be reached), on symmetry, and also on the uniformity of collision parameters (see Figs. 6 and 7). Our analytical and numerical results coincide. The distribution of waiting times has two main features. The first is a power-law decay Eq. (1), and the second are its non-analytical features (see Figs. 8 and 11). The former is due to the infinite horizon and the latter is obviously related to the periodic array of scattering centers. For example in the infinite horizon Lorentz gas, after a collision event the traveling time of the tracer particle to the next target must be larger or equal to the distance between two adjutant scattering centers (in units of unity velocity). This and other dead times implies the observed non-analytical behavior of the traveling times’ distribution. Note that for the channel model arbitrarily short times are possible, due to a “gallery whispering paths" effect, see Fig. 13 (right cell). These exhibit a very correlated motion, as a short time interval between collision events is most likely to follow by another short interval of same size, see Fig 12 (b).

This leads to the second theme of our work: can one use the distribution of time intervals between collision events to predict the time dependent probability to find particles at a given position? The answer is system dependent. The technique to make this step is based on the Lévy walk model, in dimension one for the channel and two for the Lorentz gas. This model makes the assumption of renewal, and we found it works well for the Lorentz gas and not at all for the channel. A tool to check the renewal hypothesis, from an analysis of the paths, is a correlation plot of consecutive waiting times, see Fig. 12 (b) and (d), which clearly points out the strong correlations for the channel model. However, even the failure of renewal theory does not imply the complete breakdown of the Lévy walk scheme. In this case we introduced effective (or re-normalized) parameters in the Lévy walk scheme obtained by fitting, yielding predictions that are still very useful. In fact there are general trends in the position’s distribution that are universal, and nicely predicted by the Lévy walk. These include the fat tail of the spreading packet, sharp cutoffs of the density at $|\boldsymbol{r}|=Vt$, the Kummer corrections to the Gaussian (which are certainly not small on any reasonable time scale), and the Lambert scaling. The latter is very important since it allowed us to compare finite time simulations with our theory, while the asymptotic Gaussian form (which exist for the Lorentz gas) is not witnessed due to super-slow convergence problems (see Figs. 9 and 10). One way to understand this behavior is to realize that the billiard systems are operating at a transition point between Lévy and Gauss statistics. Because of the exponent $-3$ in Eq. (1), the system is essentially behaving as if is critical in the sense of very sluggish convergence. Roughly speaking and for finite times, the packet of particles’ tails exhibit Lévy behavior (a power-law with cutoff), while the center part is Gaussian. Already in the problem of IID RVs summation, section II, we encountered a critical slowing down at this borderline case, solved by departing from the $\sqrt{N\ln(N)}$ scaling and replacing it with the Lambert approach.

To map the problem to a Lévy walk one needs to model the distribution of velocities $F_{d}(v)$. For the channel model this is rather easy, as the model is one dimensional and from symmetry we use a velocity which is either $+V$ or $-V$ with equal probability. For the Lorentz gas, a more careful analysis is needed. As we decrease the size of scatterers, we open more infinite corridors of motion. At first we have four open horizons and this leads to a cross-like shape of the spreading packet, see Fig. 2. Here, the velocity distribution in the Lévy walk scheme has a simple structure as the four directions are clearly identical from symmetry. However, when $R$ is made slightly smaller than $1/\sqrt{8}$, we open a new channel but only slightly, meaning that the effective velocity in these directions is statistically reduced compared to the original four corridors (note that we refer to the distribution of velocities in the Lévy walk, and not to the microscopic velocities of the Lorenz gas). The resulting effect is the creation of a British flag-like type of structure for the packet’s distribution, see Fig. 3. Thus, we observe a transition from one geometry to another as we vary $R$. Theoretically, this means that we assign different statistical weights to each group of horizons as in Eq. (III), which is made possible with geometrical considerations leading to the important parameter $q$ in Eq. (47).

Acknowledgements.

LZ would like to thank Jakub Ślęzak for helpful conversations. EB thanks Itzhak Fouxon for valuable discussions on the Lévy walk model. After this paper was completed a related work on the Lévy walk and Lorentz gas was presented in Itzhak . The financial support of Israel Science Foundation’s grant $1898/17$ is acknowledged (LZ and EB).

Appendix A Leading behavior of $\tilde{f}(k\to 0)$

We assume that $f(\chi)$ possesses the asymptotic behavior

[TABLE]

where $\nu>0$ and $\chi_{0}>0$ are real numbers. The Fourier transform of $f(\chi)$ is defined as

[TABLE]

Due to the evenness of $f(\chi)$, $\tilde{f}(k)$ is even and real. Therefore, throughout the following sections we assume that $k\to 0^{+}$, and use parity to find $\tilde{f}(k\to 0^{-})$. As we show below, given a positive integer $\bar{n}$, Eqs. (61-62) lead to

[TABLE]

where $C_{\bar{n}}\left[f(\chi)\right]$ is defined in the corresponding section. Using Eqs. (61) and (62), we extract all of the converging moments out of the Fourier transform integral

[TABLE]

and consider each of the cases in Eq. (63) separately.

A.1 A non-integer $\nu$

Let us assume that $\bar{n}-1<\nu<\bar{n}$. In order to find the leading behavior of the second term of Eq. (64), we consider the limit

[TABLE]

Using L’Hospital’s rule $\bar{n}$ times yields

[TABLE]

where $\Gamma(\cdots)$ is the Gamma function. We use an $\eta=k\chi$ variable change

[TABLE]

We now switch the order of limit and integration while using the asymptotic behavior (61)

[TABLE]

which proves the top row of Eq. (63).

A.2 An integer $\nu$

Let us assume that $\nu=\bar{n}$, where $\bar{n}$ is even. In order to find the leading behavior of the second term of Eq. (64), we consider the limit

[TABLE]

Using L’Hospital’s rule $\bar{n}+1$ times produces

[TABLE]

Changing the integration variable to $\eta=k\chi$ gives

[TABLE]

After switching the order of limit and integration, the integral exists as a Cauchy principal value, and we find

[TABLE]

If $\bar{n}$ is odd, we return to Eq. (68), and take the limit of $\nu\to\bar{n}$, where $\bar{n}$ is odd. We obtain

[TABLE]

To compute the next order correction for the case of an even $\bar{n}$, we calculate the following limit

[TABLE]

Using L’Hospital’s rule $\bar{n}$ times results with

[TABLE]

where $H_{\bar{n}}=\sum_{n=1}^{\bar{n}}\frac{1}{n}$ is the $\bar{n}$th harmonic number. Since $\bar{n}$ is even, the integral in Eq. (75) can be adjusted to the domain $[0,\infty)$, with $\exp(ik\chi)\to\cos(k\chi)$. We split the adjusted integral at $\chi=\chi_{0}$:

[TABLE]

We now add and subtract a $\chi_{0}^{\bar{n}}/\chi^{\bar{n}+1}$ term from $f(\chi)$ in the second integral in Eq. (76):

[TABLE]

Note that due to the asymptotics (61), the middle integral is finite when $k\to 0^{+}$. Eq. (77)’s right integral can be computed explicitly, after which the limit can be evaluated. Finally, we find

[TABLE]

where $\gamma\approx 0.5772$ is Euler’s constant. After some algebra we obtain

[TABLE]

where

[TABLE]

is a finite constant, which provides us the middle and bottom rows of Eq. (63). Plugging $\bar{n}=2$ into Eq. (80) results with Eq. (5).

Appendix B Leading behavior of $\hat{\psi}(u\to 0)$

We assume that $\psi(\tau)$ possesses the asymptotic behavior

[TABLE]

where $\nu>0$ and $\tau_{0}>0$ are real numbers. The Laplace transform of $\psi(\tau)$ is defined as

[TABLE]

As we show below, given a positive integer $\bar{n}$, Eqs. (81-82) lead to

[TABLE]

where $C_{\bar{n}}\left[\psi(\tau)\right]$ is defined in the corresponding section. Using Eqs. (81) and (82), we extract all of the converging moments out of the Laplace transform integral

[TABLE]

and consider each of the cases in Eq. (83) separately.

B.1 A non-integer $\nu$

Let us assume that $\tilde{n}-1<\nu<\tilde{n}$. In order to obtain the leading behavior of the second term of Eq. (84), we consider the limit

[TABLE]

Using L’Hospital’s rule $\tilde{n}$ times yields

[TABLE]

Changing the integration variable to $\eta=u\tau$ produces

[TABLE]

We switch the order of limit and integration, which together with the asymptotic behavior (81) gives

[TABLE]

Which proves Eq. (83)’s top row.

B.2 An integer $\nu$

Let us assume that $\nu=\bar{n}$, where $\bar{n}$ can be even or odd. In order to find the leading behavior of Eq. (84), we consider the following limit

[TABLE]

Using L’Hospital’s rule $\tilde{n}+1$, we get

[TABLE]

We change the integration variable to $\eta=u\tau$

[TABLE]

and switching the order of limit and integration, while applying the asymptotic behavior (81)

[TABLE]

To obtain the next order correction, we evaluate the following limit

[TABLE]

We use L’Hospital’s rule $\bar{n}$ times to get

[TABLE]

where $H_{\bar{n}}=\sum_{n=1}^{\bar{n}}\frac{1}{n}$ is the $\bar{n}$th harmonic number. Splitting the integral at $\tau=\tau_{0}$ yields

[TABLE]

We now add and subtract the term $\tau_{0}^{\bar{n}}/\tau^{\bar{n}+1}$ from $\psi(\tau)$ in the second integral in Eq. (95)

[TABLE]

Note that due to the asymptotics (81), the middle integral is finite when $u\to 0$. Eq. (B.2)’s right integral can be computed explicitly, after which the limit can be evaluated. Finally, we find

[TABLE]

where

[TABLE]

is a finite constant, which concludes the bottom row of Eq. (83). Plugging $\bar{n}=2$ into Eq. (97) results with Eq. (18).

Appendix C Additional steps in the derivation of Lambert scaling for the Lévy walk model

Here we portray the additional steps which were omitted in section III. Starting from Eq. (17), we expand $\hat{\psi}(u-i\boldsymbol{k}\cdot\boldsymbol{v})$ for small arguments

[TABLE]

and apply this expansion to the Montroll-Weiss Eq. (14)

[TABLE]

We use the identity

[TABLE]

and the assumed symmetry of $F(\boldsymbol{v})$ in order to simplify Eq. (99)

[TABLE]

Using $\ln(1+\epsilon^{2})\simeq\epsilon^{2}$ and $\tan^{-1}(1/\epsilon)\simeq(\pi/2)\text{sign}(\epsilon)-\epsilon$ for $\epsilon\to 0$, we discard irrelevant terms with respect to $\epsilon\sim u/kv$, so we have

[TABLE]

One may argue that this expansion breaks down when $v=0$ or alternatively when $\boldsymbol{k}\cdot\boldsymbol{v}=0$, but actually there is no problem. The former case is ruled out since our physical models have a positive constant for the speed $v=V=1$, and as such we demand that $F_{d}(v\neq 1)=0$. The latter case is ruled out since $F_{d}(\boldsymbol{v})$ covers all velocity directions of the $d$-dimensional space, by construction. Therefore, it will always contain a part parallel to $\boldsymbol{k}$, regardless of $\boldsymbol{k}$’s direction. Returning to the time domain, we obtain

[TABLE]

Substituting $\boldsymbol{\kappa}=\boldsymbol{k}\sqrt{\tau_{0}^{2}\langle v^{2}\rangle t\Omega_{d}(t)/4d\left<\tau\right>}$, where $\Omega_{d}(t)$ is a scaling function, leads to

[TABLE]

We determine the slowly varying scaling function $\Omega_{d}(t)$ by demanding that $\ln[t\Omega_{d}(t)/(4dC_{\psi}^{2}\langle\tau\rangle)]=\Omega_{d}(t)$, obtaining Eq. (III). Thus, Eq. (104) becomes

[TABLE]

where we expanded the exponential term in the second row due to the assumption of large $t$, which leads to large $\Omega_{d}(t)$. Equations (20) and (III) then follow from Eqs. (19) and (III), respectively, when combined with Eq. (C).

Appendix D Derivation of the integration boundaries

D.1 Lorentz gas model

We denote as $(x_{0},y_{0})$ the starting point on the $(0,0)$ scatterer from which we assume the particle has originated. The pair $\{b,\beta\}$ and the trio $\{x_{0},y_{0},\beta\}$ are related by a simple transformation. To obtain it, we define the two vectors

[TABLE]

where $\boldsymbol{B}$ can be seen in Fig. 6 (the red arrow). Solving the equation $\boldsymbol{B}\cdot(\boldsymbol{R}-\boldsymbol{B})=0$ and $x_{0}(0,0,R)=R$, we get the following expressions

[TABLE]

By solving

[TABLE]

together with $\tau_{1,0}^{*}(0,0,R)=1-2R$, we obtain Eqs. (IV.1) and (29). In order to extract the integration boundaries (IBs) of $\beta$ and $b$ out of the inequalities Eqs. (30-IV.1), we notice that there are two classes of trajectories which reach the $(n,1)$ scatterer, where $n\geq 0$. The first class’ IBs are dictated by the origin and target scatterers. The second class’ IBs are governed by the $(1,0)$ and $(n-1,1)$ scatterers [this class does not exist for the $(0,1)$ circle]. We split the $\beta$ domain into two parts marked (I) and (II), each corresponds to a class of trajectories, and denote the separator angle between them with $\beta^{\rm sep}_{n,1}(R)$. Each of $\beta$’s subdomains is associated with a different expression for $b$’s IBs, which we denote (i) and (ii). However, when $\beta=\beta^{\rm sep}_{n,1}(R)$ the two expressions coincide. Therefore, $\beta^{\rm sep}_{n,1}(R)$ can be found by comparing the lower (or upper) IB of (i) to that of (ii). Finding the top/bottom IBs of $\beta$ relies on a similar idea, namely using the expressions for $b$’s IBs. When $\beta=\beta^{\rm min}_{n,1}(R)$, the associated $b$ subdomain (i) shrinks to zero. Thus, $\beta^{\rm min}_{n,1}(R)$ can be found by setting $b^{\rm min}_{n,1}(\beta,R)=b^{\rm max}_{n,1}(\beta,R)$ on the (i) expression. Identically, $\beta^{\rm max}_{n,1}(R)$ is found by setting $b^{\rm min}_{n,1}(\beta,R)=b^{\rm max}_{n,1}(\beta,R)$ on the (ii) expression. Starting with $n=0$, the upper IB of $\beta$ is set to $\pi/2$, which is possible due to symmetry. This is done in order to balance out the excess of distant neighbors groups over nearest and next to nearest neighbors groups across the lattice, which is demonstrated in Fig. 7. Writing Eq. (30) in its explicit form, we get

[TABLE]

Since $\beta\leq\pi/2$ and $|b|\leq R$, the upper IB in Eq. (109) is set to $R$. For $\beta$’s lower IB, we equate $b$’s IBs, $\cos(\beta)-R=R$, getting $\cos^{-1}(2R)$. Thus, we obtain for $n=0$

[TABLE]

We continue with $n=1$. Again, due to symmetry we can truncate $\beta$’s upper IB to $\pi/4$. Equation (IV.1) in its explicit form combined with $|b|\leq R$ then yields

[TABLE]

Here, the $\beta$ domain is sectioned into two parts, as said before: (I) for which the (ii) inequality in Eq. (111) is trivially fulfilled, and (II) in which it needs to be upheld. Equating the IBs of (i) in Eq. (111) yields the lower IB of $\beta$. We obtain

[TABLE]

The separator angle between (I) and (II) can be found by comparing the lower (or upper) $b$ IB of (i) to that of (ii), yielding $\beta=\cos^{-1}(2R)$. Therefore we find for $n=1$

[TABLE]

and

[TABLE]

Finally, Eq. (IV.1) in its explicit form supply

[TABLE]

where we already implemented $|b|\leq R$ to (i)’s top IB. This time, there is no need to set the upper IB of $\beta$ manually. For the bottom IB of $\beta$’s subdomain (I), one equates Eq. (115)’s (i) IBs, namely $\cos(\beta)-n\sin(\beta)=2R$. By dividing this equality with $\sqrt{n^{2}+1}$ and using basic trigonometry, we have a general form of the lower IB of $\beta$, Eq. (33). The same is done for the upper IB of $\beta$’s subdomain (II), where (ii) of Eq. (115) is used, yielding $\cos(\beta)-(n-2)\sin(\beta)=2R$, which results with Eq. (34). The separator angle can be found by equating the top/bottom IBs of (i) to (ii), producing $\cos(\beta)-(n-1)\sin(\beta)=2R$ and Eq. (36). We finish this subsection with the case of $R=0.3$, namely four open corridors. As said in section IV, the scatterer $(2,1)$ is now shared by the horizontal and diagonal directions, thus we adjust its upper $\beta$ IB. When $\beta$ achieves its maximal value for the horizontal direction stripe, $b$’s lower IB must be equal to $-R$. Thus we set $b^{\rm min}_{2,1}(\beta,R)=-R$ in Eq. (35), and obtain Eq. (39). Similarly to the horizontal case, there are two classes of trajectories to reach the $(m+1,m)$ scatterer. The first class’ IBs are governed by the origin and target scatterers, and the second class’ IBs are dictated by the $(1,1)$ and $(m,m-1)$ circles [this class does not exist for the $(2,1)$ target scatterer]. Again, we split $\beta$’s domain into two subdomains (I) and (II), however this time they switch places, i.e. the (I) part is of higher $\beta$ values than the (II) part. Equation (40) in its explicit form reads

[TABLE]

where we set the lower $b$ IB to $-R$ since $\beta\geq\sin^{-1}(2R)$ and $|b|\leq R$. For the upper IB of $\beta$ we equate the IBs of $b$ in Eq. (116) to each other, namely $-R=\cos(\beta)-2\sin(\beta)+R$. Thus we obtain for $m=1$

[TABLE]

For $m>1$, we write Eq. (IV.1) in its explicit form

[TABLE]

where we set the bottom IB of (i) to $-R$, as before. The upper $\beta$ IB is found by equating (i)’s top IB to $-R$ in Eq. (118), and the bottom $\beta$ IB is found by equating (ii)’s top and bottom IBs to each other. The separator angle is found in the same manner as for the horizontal case. Dividing the resulted expressions with $\sqrt{m^{2}+(m+1)^{2}}$ and using basic trigonometry, we obtain Eqs. (IV.1-45).

D.2 Stadium channel model

We denote as $(x_{0},y_{0})$ the starting point on the $(0,0)$ stadium from which we assume the particle has originated. The pair $\{a,\alpha\}$ and the trio $\{x_{0},y_{0},\alpha\}$ are related by a simple transformation. To obtain it, we define the two vectors

[TABLE]

where $\boldsymbol{A}$ can be seen in Fig. 6 (the red arrow). Solving the equation $\boldsymbol{A}\cdot(\boldsymbol{R}-\boldsymbol{A})=0$ and $y_{0}(0,0)=-1$, we get the following expressions

[TABLE]

By solving

[TABLE]

together with $\tau_{0,0}^{*}(0,0)=2$, we obtain Eqs. (IV.2) and (IV.2). Simplifying Eq. (51) and using $-\pi/2\leq\alpha\leq\pi/2$, we get Eq. (IV.2). In order to extract $\alpha$’s lower IB out of the inequalities Eq. (52), we notice that as with the Lorentz gas model, when $\alpha$ hits its bottom IB, the $a$ domain vanishes. Simplifying Eq. (52), we have

[TABLE]

and thus we have $D\cos(\alpha)-\sin(\alpha)=\sin(\alpha)$ for $\alpha$’s lower IB, which solves to $\alpha_{0,D}^{\rm min}=\tan^{-1}(D/2)$. Finally, we notice that there are two classes of trajectories which reach the target stadiums of $n>0$, as with the Lorentz gas case. The first class’ IBs are dictated by the origin and target semicircles. The second class’ IBs are governed by the origin and $(2n-2,D)$ stadium walls [this class does not exist for the trajectories ending with the origin or $(0,D)$ semicircles]. We split the $\alpha$ domain into two parts marked (I) and (II), each corresponds to a class of trajectories, and denote the separator angle between them with $\alpha^{\rm sep}_{2n,D}$. Each of $\alpha$’s subdomains is associated with a different expression for $a$’s IBs, which we denote (i) and (ii). However, when $\alpha=\alpha^{\rm sep}_{2n,D}$ the two expressions coincide. Therefore, $\alpha^{\rm sep}_{2n,D}$ can be found by comparing the lower (or upper) IB of (i) to that of (ii). Finding the top/bottom IBs of $\alpha$ relies on a similar idea, namely using the expressions for $a$’s IBs. When $\alpha=\alpha^{\rm min}_{2n,D}$, the associated $a$ subdomain (i) shrinks to zero. Thus, $\alpha^{\rm min}_{2n,D}$ can be found by setting $a^{\rm min}_{2n,D}(\alpha)=a^{\rm max}_{2n,D}(\alpha)$ on the (i) expression. Identically, $\alpha^{\rm max}_{2n,D}$ is found by setting $a^{\rm min}_{2n,D}(\alpha)=a^{\rm max}_{2n,D}(\alpha)$ on the (ii) expression. Simplifying Eq. (IV.2)’s top row, we get

[TABLE]

This class of trajectories has its upper/lower $a$ IB dominated by the origin/target stadium, and thus the top/bottom IB for this class is taken from Eq. (123)’s first/second inequality, such that

[TABLE]

For the inequalities of Eq. (IV.2)’s bottom row, we get after simplification

[TABLE]

This class of trajectories has its lower/upper $a$ IB dominated by the origin/$(2n-2,D)$ stadium, and thus the bottom/top IB for this class is taken from Eq. (123)’s first/second inequality, such that

[TABLE]

Using the IBs of $a$ in (i) and (ii) to extract the IBs of $\alpha$ in the way described above, we obtain Eqs. (55-58).

Appendix E Calculations of $\tau_{0}$, $\langle\tau\rangle$, and $C_{\psi}$

E.1 Lorentz gas model

Using Eq. (28) together with the integration boundaries Eqs. (33-36), one can calculate $\psi(\tau)$ using a computational program like Mathematica and extract the constants $\tau_{0}$, $\langle\tau\rangle$, and $C_{\psi}$ out of it. However, we found that an analytical expression for $\tau_{0}$ can be calculated. It follows from its definition Eq. (1) that

[TABLE]

Plugging Eq. (28) with two infinite corridors (i.e. $q=0$) into Eq. (127), we get

[TABLE]

where $\text{H}(\cdots)$ is the Heaviside step function, and the factor of $8$ arise from symmetry, see Fig. (7). Since $\tau_{n,m}^{*}(b,\beta,R)$ is the traveling distance to the $(n,m)$ scatterer, it obeys $\tau_{n,1}^{*}(b,\beta,R)\simeq\sqrt{n^{2}+(1-2R)^{2}}\simeq n$ when $n$ is large, thus the Heaviside function truncates the sum in Eq. (128) at $n=T$, and we have

[TABLE]

Equation (127) suggests that the integral over $\tau^{3}\psi(\tau)$ grows linearly with $T$, and as such we expect that the sum in Eq. (129) will behave similarly with $T$. We therefore write

[TABLE]

Now the integrals can be easily performed. After evaluating the limit we get a closed expression for $\tau_{0}$, Eq. (37). The remaining constants $\langle\tau\rangle$ and $C_{\psi}$ are defined via integrals over $\psi(\tau)$ rather than by a limit operation, and as such they cannot be obtained using end terms as we just did. Even though, one is not required to calculate $\psi(\tau)$, but can use a simpler tactic. To calculate the mean time between collisions, we use its definition and plug inside Eq. (28)

[TABLE]

To achieve a designated precision, one can simply truncate the sum at a large enough $T$. For $T=500$ we obtain $\langle\tau\rangle\approx 0.62155$. For $C_{\psi}$ we find the following formula out of Eq. (18)

[TABLE]

Plugging Eq. (28) into the above expression and allowing the Heaviside function to truncate the sum at $n=T$ yields

[TABLE]

This equation converges rather slowly due to the logarithmic term, hence we need to accelerate its convergence rate. To do that, we use the identity

[TABLE]

and write, taking the $n=0$ summand out of the sum

[TABLE]

We approximate $\tau_{n,1}^{*}(b,\beta,R)$ for large $n$ as before and see, after performing the integrals, that the summands behave for $n\gg 1$ as

[TABLE]

where $B_{1},B_{2},B_{3}$ are some $R$-dependent coefficients. This behavior suggests that the partial sum $S_{T}$ goes like

[TABLE]

for large $T$, where $l$ is the desired limit. This implies that the Richardson extrapolation method can be used, see Eq. (8.1.16) in Ref. Bender . For $T=40$, we get $C_{\psi}\approx 4.4802\times 10^{-4}$ by summing terms up to $n=T$ and extrapolating over the $S_{37},S_{38},S_{39},S_{40}$ partial sums.

Extracting the needed constants for four open horizons, namely $q>0$, is done in a similar fashion, but this time we have $1/\sqrt{20}\leq R<1/\sqrt{8}$. Plugging Eq. (28) into Eq. (127) gives

[TABLE]

Using geometrical considerations, this time we have $\tau^{*}_{m+1,m}(b,\beta,R)\simeq[2(m+1/2)^{2}+(1/\sqrt{2}-2R)^{2}]^{1/2}\simeq\sqrt{2}m$ for large $m$, so Eq (138)’s sum is truncated at (the closest integer to) $T/\sqrt{2}$

[TABLE]

After loosing the sum as was done for the two open horizons case, Eq. (139) becomes

[TABLE]

which then yields Eq. (46). We are left with calculating the remaining parameters for the flag case. For the mean time between collisions we have

[TABLE]

Notice that if one truncates the first sum in Eq. (141) at $T$, one then needs to truncate the second sum at $T/\sqrt{2}$. Using $T=500$, we obtain $\langle\tau\rangle\approx 1.1947$. For $C_{\psi}$ we write

[TABLE]

and consequently

[TABLE]

The behavior displayed in Eq. (137) was checked to be valid for both sums in Eq. (E.1). Thus, we take $T=40$ and use Richardson extrapolation for each sum separately by extrapolating over the last four partial sums (for the second sum we use terms for which $m\leq 28\approx T/\sqrt{2}$), and obtain $C_{\psi}\approx 1.5250\times 10^{-2}$.

E.2 Stadium channel model

We use an identical way as for the Lorentz gas to compute $\tau_{0}$ analytically. Since $\tau_{2n,D}^{*}(a,\alpha)$ is the traveling distance to the $n$th upper semicircle, it is clear that $\tau_{2n,D}^{*}(a,\alpha)\simeq\sqrt{(2n+2)^{2}+D^{2}}$ when $n$ is large. Plugging Eq. (49) into Eq. (127) and employing similar manipulations as before, we have

[TABLE]

Since $\sqrt{(2n+2)^{2}+D^{2}}\simeq 2n$ for large $n$, the Heaviside function truncates the sum in Eq. (144) at $n=T/2$. Thus

[TABLE]

and consequently

[TABLE]

Now the integrals can be easily performed. After evaluating the limit we get a closed expression for $\tau_{0}$, Eq. (59). For the mean time between collisions, we write

[TABLE]

We truncate the sum at $T=500$, and obtain $\langle\tau\rangle\approx 2.57016$. For $C_{\psi}$ we have

[TABLE]

Substituting $\tau_{2n,D}^{*}(a,\alpha)\simeq\sqrt{(2n+2)^{2}+D^{2}}$ for $n\gg 1$ and evaluating the integrals, we observe the following behavior

[TABLE]

where $A_{1},A_{2},A_{3}$ are some $D$-dependent coefficients. This behavior suggests that the partial sum $S_{T}$ goes like

[TABLE]

for large $T$, where $l$ is the desired limit. Once again, we employ the Richardson extrapolation method. For $T=60$, we get $C_{\psi}\approx 1.0903\times 10^{-5}$ by summing terms up to $n=T/2$ and extrapolating over the $S_{27},S_{28},S_{29},S_{30}$ partial sums.