On the free path length distribution for linear motion in an n-dimensional box
Samuel Holmin Pär Kurlberg Daniel Månsson


Trace of the Twisted Heisenberg Category
    Michael Reeks
    December 30, 2023
========================================




  We consider the distribution of free path lengths, or the distance
  between consecutive bounces of random particles, in an
  n-dimensional rectangular box.  If each particle travels a distance
  R, then, as R →∞ the free path lengths coincides with
  the distribution of the length of the intersection of a random line
  with the box (for a natural ensemble of random lines) and we
  give an explicit formula
  (piecewise real analytic)
  for the probability density function in dimension two and three.

  In dimension two we also consider a closely related model where each
  particle is allowed to bounce N times, as N →∞, and give
  an explicit (again piecewise real analytic)
  formula for its probability density function.

  Further, in both models we can recover the side lengths of the box
  from the location of the discontinuities of the probability density
  functions.










§ INTRODUCTION

We consider billiard dynamics on a rectangular domain, i.e., point
shaped “balls” moving with linear motion with specular reflections at
the boundary, and similarly for rectangular box shaped domains in
three dimensions.  We wish to determine the distribution of free path
lengths of ensembles of trajectories defined by selecting a starting
point and direction at random.


The question seems quite natural and interesting on its own, but we
mention that it originated from the study of electromagnetic fields in
“reverberation chambers” under the assumption of highly directional
antennas <cit.>.  Briefly, the connection is as
follows (we refer to the forthcoming paper <cit.> for
more details): given an ideal highly directional antenna and a highly
transient signal,
then the wave pulse dynamics is essentially the same as a point shaped billiard
ball traveling inside a chamber, with specular reflection at the
boundary.  Signal loss is dominated by (linear) “spreading” of the
electromagnetic field and by absorption occurring at each interaction
(“bounce”) with the walls.
The first simple model we use in this paper neglects absorption
effects, and models signal loss from spreading by simply terminating
the motion of the ball after it has travelled a certain large
distance.  The second model only takes into account signal loss from
absorption, and completely neglects spreading; here the motion is
terminated after the ball has bounced a certain number of times.

We remark that the distribution of free path lengths is very well
studied in the context of the Lorentz gas — here a point particle
interacts with hard spherical obstacles, either placed randomly, or
regularly on Euclidean lattices; recently quasicrystal configurations
have also been studied (cf. <cit.>.)




Let R>0 be large and let a rectangular n-dimensional box
K⊆^n be given, where n≥ 2. We send off a large number M>0 of particles,
each with a random initial position p^(i)∈ K chosen with respect
to a given probability measure μ on K, and each with a uniformly
random initial direction v^(i)∈𝕊^n-1={x∈^n: 
x=1}, i=1,…,M, for a total distance R each. 
Each particle travels along straight lines, changing direction
precisely when it hits the boundary of the box, where it reflects
specularly. We record the distance 
travelled between each pair of consecutive bounces for each particle.
(Note in particular that we obtain more bounce lengths from some
particles than from others.) 


Let X_M,R be the uniformly distributed random variable on
this finite set of bounce lengths of all the particles.
More precisely, a random sample of X_M,R is obtained as follows:
first take a random i.i.d.  sample of points (with respect to the
measure μ) p^(1), …, p^(M)∈ K, and a random sample of
directions v^(1),…,v^(M)∈𝕊^n-1 (with respect to
the uniform measure).  Each pair (p^(i),v^(i)) then defines a
trajectory T^i of length R, and each such trajectory gives rise
to a finite multiset B^i of lengths between consecutive
bounces.  Finally, with B = _i=1^M B^i denoting the
(multiset) union of bounce length multisets B^1,…,B^M, we
select an element of B with the uniform distribution.  (That is, with
1_B denoting the integer valued set indicator function for B,
and B' = { x : 1_B(x) ≥ 1 } we select the element b ∈ B'
with probability 1_B(b)/∑_x ∈ B' 1_B(x).)


We are interested in the distribution of X_M,R for large M and
R, and this turns out to be closely related to a model arising from
integral geometry.
Namely, let ℓ denote the unique (up to a constant) translation-
and rotation-invariant measure on the set of directed lines ℓ in
^n, and consider the restriction of this measure to the set of
directed lines ℓ intersecting K, normalized such that it
becomes a probability measure. Denote by X the random variable
Xlength(ℓ∩ K) where ℓ is chosen at
random using this measure. 


For any dimension n≥ 2, and for any distribution μ on the
starting points, the random variable X_M,R converges in distribution
to the random variable X, as we take R→∞ followed by taking
M→∞, or vice versa.


The mean free path length has a quite simple geometric interpretation. We have
  
    X = 2π|𝕊^n-1||𝕊^n|Vol(K)Area(K)
         = 2√(π)·Γ(n+1/2)Γ(n/2)Vol(K)Area(K)

  where Area(K) is the (n-1)-dimensional surface
  area of the box K, Vol(K) is the volume of the
  box K, Γ is the gamma function, and where
  |𝕊^n-1|=2π^n/2/Γ(n/2) is the (n-1)-dimensional surface area
  of the sphere 𝕊^n-1⊆^n.
  The formula in (<ref>) has been proven in a more general setting earlier (see formula (2.4) in <cit.>); for further details, see Section <ref>. For the convenience of the reader we give a short proof of formula (<ref>) in our setting in Section <ref>.


Throughout the paper, we will write pdf_Z and cdf_Z for
the probability density function and the cumulative distribution
function of Z, respectively, for random variables Z. 
We next give explicit formulas for the probability density function of X
in dimensions two and three.



  For a box of dimension n=2 with side-lengths ≤, the probability density function of X is given by
  
    pdf_X(t) = 1/+·
          1,   if t<,
    ^2t^2√(t^2-^2),   if <t<
    
          -1 + 1t^2( ^2√(t^2-^2) + ^2√(t^2-^2)),   if ,<t.

  for 0<t<√(^2+^2).




  We note that the probability density function in Theorem <ref>
  is analytic on all open subintervals of (0, √(^2+^2)) not
  containing  or . Moreover, it is constant on the interval
  (0,min(,)) and has singularities of type (t-)^-1/2 and
  (t-)^-1/2 just to the right of  and , respectively.
  See Figure <ref> for more details.  For an
  explanation of these singularities, see
  Remark <ref>.



 






  For a box of dimension n=3 with side-lengths ,,, the probability density function of X is given by
  
    pdf_X(t) = F(,,,t) + F(,,,t) + F(,,,t)3 π  t^3 (++)

  where F is the piecewise-defined function given by
  
    F(,,,t) = t^3(8-3t)

  for 0<t<, and by
  
    F(,,,t) =
              (6t^4-^4+6 π^2 )-4(+)√(t^2-^2)(^2+2t^2)

  for <t<√(^2+^2),
  and by
  
    F(,,,t) = 6 π^2  +^4-3 t^4-6 ^2 ^2 + 
    √(t^2-^2-^2) 4( ^2 + ^2 +2 t^2)+
    
                +4√(t^2-^2)(^2+2t^2)-12 ^2 ·arctan(√(t^2-^2-^2))+
    
                -4√(t^2-^2)(^2+2t^2)-12 ^2 ·arctan(√( t^2-^2-^2)/)

  for √(^2+^2)<t<√(^2+^2+^2).







  We note that the probability density function in Theorem <ref> is analytic on all open subintervals of (0, √(^2+^2+^2)) not containing any of the points
  
    ,,,√(^2+^2), √(^2+^2), √(^2+^2).

  Moreover, it is linear on the interval (0,min(,,)) and has
  positive jump discontinuities at the points ,,. At the
  points {√(^2+^2), √(^2+^2),
  √(^2+^2)}∖{,,}, it is continuous and differentiable.


Note that the probability distribution X_M,R gives a larger
“weight” to some particles than others, since some particles get
more bounces than others for the same distance R. One could also
consider a similar problem where we send off each particle for a
certain number N>0 of bounces, and then consider the limit as
M→∞ followed by taking the limit N→∞, where M is
the number of particles. This would give each particle the same
“weight”.  Denote the finite version of this distribution by
Y_M,N and its limit distribution as M→∞ and then N→∞ by Y.  With
regard to the previous discussion about signal loss, we call the limit
distribution X of X_M,R the spreading model and we call
the limit distribution of Y_M,N the absorption
  model. Determining the probability density function of the absorption model appears
to be the more difficult problem, and we give a formula only in dimension
two:



  For a box of dimension n=2 with side-lengths ≤, the random variable Y_M,N converges in distribution
  to the random variable Y, as we take M→∞ followed by taking N→∞, where
  the probability density function pdf_Y(t) is given by
    
    2/π(
          2 (+)/(^2+^2)-2 /(^2+^2)^3/2(tanh^-1(/√(^2+^2)) + tanh^-1(/√(^2+^2)))
          )

    for 0<t<,, and by
    
    2/π(
          (- √(t^2-^2)) /t  (+√(t^2-^2)) √(t^2-^2) +
          2  +2 t - 2√(t^2-^2)/t (^2+^2)+
    2 (-tanh
           ^-1(t/√(^2+^2))+tanh ^-1(√(t^2-^2)√(^2+^2)/t ) - tanh^-1(/√(^2+^2)))/(^2+^2)^3/2)

    for <t<, and by
    
    2/π(
          (- √(t^2-^2))/t  (+√(t^2-^2)) √(t^2-^2) 
          +
          (- √(t^2-^2))/t  (+√(t^2-^2)) √(t^2-^2) 
          +
          22- √(t^2-^2) - √(t^2-^2)/t (^2+^2)+ 
    2 (-2tanh
           ^-1(t/√(^2+^2))+tanh ^-1(√(t^2-^2)√(^2+^2)/t ) +
            tanh ^-1(√(t^2-^2)√(^2+^2)/t )
            )/(^2+^2)^3/2)

    for ,<t<√(^2+^2).




See Figure <ref> for a comparison between the probability density functions for the two different models in dimension 2.

  It is not a priori obvious that the two limit distributions should
  differ, and it is natural to ask how much, if at all, they
  differ. We start by remarking that the expression for
  pdf_Y(t) does not simplify into the expression for
  pdf_X(t); indeed, for (,)=(1,2) we have
  pdf_X(t)=1/3 but pdf_Y(t)≈ 0.32553 on the
  interval (0,1).  For very skew boxes, with =1 and
  →∞, it is straightforward to show that

    pdf_Y(/2)/pdf_X(/2)→∞

as →∞.




 §.§ Discussion



Given a closed convex subset C ⊂^n with nonempty interior
it is possible to define a natural
probability measure on the set of lines in ^n that have nonempty
intersection with C.  The expected length of the intersection of
a random line is then, up to a constant that only depends on n, 
given by Vol(C)/Area(C); this is known
as Santalo's formula 
in the integral geometry 
and geometric probability literature (cf. <cit.>).


A billiard flow on a manifold M with boundary ∂ M gives
rise to a billiard map (roughly speaking, the phase space Ω is
then the collection of inward facing unit vectors v at each point
x ∈∂ M).  Given (x,v) ∈Ω we define the
associated free path as the distance the billiard particle, starting
at x in the direction v, covers before colliding with ∂ M
again.  As the billiard map carries a natural probability measure
ν we can view the free path as a random variable, and the mean
free path is then just its expected value.  Remarkably, the mean free
path (again up to a constant that only depends on the dimension) is
then given by Vol(M)/Vol(∂ M) — even for non-convex
billiards.  This was deduced in the seventies at the Moscow seminar on
dynamical systems directed by Sinai and Alekseev but was never
published and hence rederived by a number of researchers.  For further
details and an interesting historical survey, see Chernov's paper
<cit.>.

In spirit our methods are closely related to the ones used by
Barra-Gaspard <cit.> in
their study of the level spacing distribution for quantum 
graphs, and this turns out to be given by the distribution of return
times to a hypersurface of section of a linear flow on a torus.

In particular, for graphs with a finite number of disconnected bonds
of incommensurable lengths, the hypersurface of section is the
“walls” of the torus, and the level spacings of the quantum graph is
exactly the same same as the free path length distribution in our
setting when all particles have the same starting velocity.  (In
particular, compare the numerator in (<ref>) for v
fixed with <cit.>.)

In <cit.>, Marklof and
Strömbergsson used the results by Barra-Gaspard to determine the gap
distribution of the sequence of fractional parts of
{log_b n }_n ∈^+.  The gap distribution depends on
whether b is trancendental, rational or algebraic;  quite
remarkably the density function P(s) for these gaps share a number
of qualitative features with the density function
pdf_X(s) for free paths in our setting.  Namely, the
density functions both have compact support and are smooth apart from
a finite number of jump discontinuities.  Further, in some cases the
density function is constant for s small; compare
Figure <ref> (here d=2) with
<cit.> (here
b=√(10)).  However, there are some important differences: for
P(s), left and right limits exist at the jump discontinuities,
whereas for d=2, the right limit of pdf_X(s) is
+∞ at the jumps (cf. Figure <ref>.)
Further, despite appearences, P(s) is not 
linear near s=0
(cf. <cit.>
corresponding to b=e) whereas for d=3, pdf_X(s)
is indeed linear near s=0 (cf. Figure <ref>).





 §.§ Acknowledgements

We would like to thank Z. Rudnick for some very helpful discussions,
especially for suggesting the connection with integral geometry.  We
also thank J. Marklof for bringing references
<cit.> to our attention.

S.H. was partially supported by a grant from the
  Swedish Research Council (621-2011-5498).

P.K. was partially supported by grants from the Göran Gustafsson
  Foundation for Research in Natural Sciences and Medicine, and the
  Swedish Research Council (621-2011-5498).



§ PROOF OF THEOREM <REF>

In this section, we prove Theorem <ref>. 
For notational simplicity, we give the proof in dimension three; the general proof for n≥ 2 dimensions is analogous.

Given a particle with initial position p and initial direction v, let N_R,p,v be the number of bounce lengths we get from that particle as it has travelled a total distance R>0, and let N_R,p,v(t) be the number of such bounce lengths of length at most t≥ 0. The uniform probability distribution on the set of bounce lengths of M particles with initial positions p^(1),…,p^(M) and initial directions v^(1),…,v^(M) has the cumulative distribution function

    cdf_X_M,R(t)= 
    ∑_i=1^M N_R,p^(i),v^(i)(t)∑_i=1^M
      N_R,p^(i),v^(i) 
    = 1/M∑_i=1^M
      N_R,p^(i),v^(i)R
      N_R,p^(i),v^(i)(t)N_R,p^(i),v^(i)1/M∑_i=1^M  N_R,p^(i),v^(i)R
    .

(Note that the denominator is uniformly bounded from below, which follows from equation (<ref>) below.)
By the strong law of large numbers, the function (<ref>) converges almost
surely to

    ∫_K∫_𝕊^2N_R,p,vR
              N_R,p,v(t)N_R,p,v S(v)μ
              p∫_K∫_𝕊^2N_R,p,vR 
              S(v)μ p

as M→∞, where dμ is the probability measure with which we
choose the starting points, and S is the surface area measure on
the sphere 𝕊^2.  By symmetry, we may restrict the inner
integrals to 𝕊^2_+{(,,)∈𝕊^2:
,,> 0}. 
We now look at the limit of (<ref>) as R→∞, and
we note that since the integrands are uniformly bounded, we may move
the limit inside the integrals by the Lebesgue dominated convergence
theorem. Fix one of the integrands, and denote it by f(R,p,v,t).
We will show that its limit g(p,v,t)lim_R→∞
f(R,p,v,t) exists for all t and all directions v∈𝕊^2. 
Moreover, if p^(i) and v^(i) denote random variables corresponding to
an initial position and an initial direction, respectively, as above, then

    h(p^(i),v^(i),t) := lim_R →∞N_R,p^(i),v^(i)R
      N_R,p^(i),v^(i)(t)N_R,p^(i),v^(i)

is a random variable with
finite variance (and similarly for the terms in the denominator of
(<ref>); in particular recall it is uniformly bounded
from below), 
and thus the strong law of large numbers gives that
the limit of (<ref>) as R→∞, and then M →∞
almost surely equals (<ref>).   This shows that 
lim_M→∞lim_R→∞cdf_X_M,R(t) exists
almost surely
and is equal to 
lim_R→∞lim_M→∞cdf_X_M,R(t).

Consider a particle with initial position p and initial direction
v=(,,)∈𝕊^2_+. By “unfolding” its motion with
specular reflections on the walls of the box to the motion along a
straight line in ^n — see Figure <ref> for a 2D
illustration — we see that the particle's set of bounce lengths is
identical to the set of path lengths between consecutive intersections
of the straight line segment {p+tv: 0≤ t≤ R} with any of the
planes x=n, y=n, z=n, n∈. Thus we see that

    N_R,p,v = R/+R/+R/+(1)

for
large R, and therefore

    N_R,p,vR→/+/+/

as R→∞.




Now project the line {p+tv:0≤ t≤ R} to the torus ^3/Λ where Λ={(n_1,n_2,n_3): n_1,n_2,n_3∈} and let us identify the torus with the box K; see Figure <ref>. Each bounce length corresponds to a line segment which starts in one of the three planes x=0, y=0 or z=0 and runs in the direction v to one of the three planes x=, y= or z=. There are
R/+(1)
line segments which start from the plane z=0, and thus the probability that a line segment starts from the plane z=0 is

    //+/+/

as R→∞.  By the ergodicity of the linear flow on tori (for
almost all directions),
the starting
points of these line segments become uniformly distributed on the
rectangle [0,]×[0,]×{0} for almost all
v∈𝕊^2_+ as R→∞; from here we will assume that
v is such a direction, and we will ignore the measure zero set of
directions for which we do not have ergodicity. 
Consider one of these line segments and denote its length by T and
its starting point by (x_0,y_0,0). For an arbitrary parameter
t≥ 0, we have T≤ t if and only if t≥-x_0 or
t≥-y_0 or t≥; the starting points
(x_0,y_0)∈[0,]×[0,] which satisfy this are precisely
those outside the rectangle [0,-t]×[0,-t] assuming
that t≤ and otherwise it is the whole rectangle
[0,]×[0,]. The area of that region is

    -(-t)(-t)

if ≥ t, ≥ t, ≥ t and otherwise it is . Since the starting points (x_0,y_0) are uniformly distributed in the rectangle [0,]×[0,] as R→∞, it follows that the probability that T≤ t is

    1-(-t)(-t)χ(≥ t, ≥ t, ≥ t),

where χ(P) is the indicator function which is 1 whenever the condition P is true, and 0 otherwise.
We get analogous expressions for the case when a line segment starts
in the plane x=0 or y=0 instead. Thus the
proportion
   of all line 
segments with length at most t as R→∞ is 

  
    lim_R→∞N_R,p,v(t)N_R,p,v   =
      	//+/+/(1-(-t)(-t)χ(≥ t, ≥ t, ≥ t))+
       =.//+/+/(1-(-t)(-t)χ(≥ t, ≥ t, ≥ t))+
       =.//+/+/(1-(-t)(-t)χ(≥ t, ≥ t, ≥ t))


which can be written

    1-χ(≥ t, ≥ t, ≥ t)
              (/+/+/)×
    ×(
      		(-t)(-t)+
      		(-t)(-t)+
      		(-t)(-t)
      	).

Recognizing that both integrands (<ref>) and (<ref>) are independent of the position p, we see that the limit of (<ref>) as R→∞ may be written as


    lim_R→∞lim_M→∞cdf_X_M,R(t) =
    	1-
    	1
    	∫_𝕊^2_+
    		(++)
    	 S(v) ×
    ×∫_
    			v∈𝕊^2_+
    ≤/t 
    ≤/t 
    ≤/t   
    		(
    		(++)
    		-2t(++)
    		+3t^2
    	)
    	 S(v)

for all t>0.
The corresponding formula in n dimensions is given by
  
    lim_R→∞lim_M→∞cdf_X_M,R(t) =
        1-
        ∫_
          v∈𝕊^n-1_+
    
          v_i ≤_i/t 
     for i=1,…,n   (
              ∑_i=1^n
                v_i ∏_j≠ i(a_i-tv_j)
            )
           S(v)
        (∏_i=1^n _i) ∫_𝕊^n-1_+( ∑_i=1^n v_i/_i)
         S(v)

for all t>0, where the side-lengths of the box K are
_1,…,_n and S is the surface area measure on 𝕊^n-1_+∩ [0,∞)^n. (The denominator can be given explicitly by
using Lemma <ref> below.) 

We have thus proved that the random variable X_M,R converges
in distribution to a random variable with probability density function
given by (<ref>) as we take M→∞ followed by
taking R→∞, or alternatively, first taking R→∞ followed by taking M→∞.
It remains to prove that this distribution agrees with
the distribution of the random variable X defined in the
introduction. 



 §.§ Integral geometry

We start by recalling some standard facts from integral geometry (cf.
<cit.>.)

The set of directed straight lines ℓ in ^3 can be
parametrized by pairs (v,q) where v∈𝕊^2 is a unit
vector pointing in the same direction as ℓ and q∈ v^ is
the unique point in ℓ which intersects the plane through the
origin which is orthogonal to v. The unique translation- and
rotation-invariant measure (up to a constant) on the set of directed
straight lines in ^3 is ℓ A(q)  S(v) where
A is the surface measure on the plane through the origin
orthogonal to v∈𝕊^2, and S is the surface area
measure on 𝕊^2. 

Consider the set L_,, of directed straight lines in ^3
which intersect the box K. Now, since ++ is
the area of the projection 
of the box K onto the plane v^ for
v∈𝕊^2_+, it follows that the total measure of L_a,b,c with
respect to ℓ is

    C_,, 8∫_𝕊^2_+(++) S(v) = 2π
            (++)

where we used symmetry, and the integral may be evaluated by switching to spherical coordinates. It follows that ℓ / C_,, is a probability measure on the set of directed lines intersecting the box L_a,b,c. Let ℓ be a random directed line with respect to this measure, and define the random variable Xlength(ℓ∩ K), as in the introduction. Let us determine the probability that X≤ t for an arbitrary parameter t≥ 0. By symmetry it suffices to consider only directed lines with v∈𝕊^2_+.
The set of all intersection points between the rectangle
[0,]×[0,]×{0} and the lines ℓ with X≤ t
and direction v∈𝕊^2_+ has area
-(-t)(-t)χ(≥ t,≥ t,≥ t), as
in (<ref>), and its projection onto the plane v^ has
area 

    [-(-t)(-t)χ(≥ t,≥ t,≥ t)].

By symmetry it follows that the area of the set of directed lines
ℓ∈ L_,, with X≤ t and direction v∈𝕊^2_+ projected down to v^ is 

    U(v,t)    [-(-t)(-t)χ(≥ t,≥ t,≥ t)] + 
       .[-(-t)(-t)χ(≥ t,≥ t,≥ t)] + 
       .[-(-t)(-t)χ(≥ t,≥ t,≥ t)],

and it follows that

    Prob[X≤ t] = 1C_,,∫_X≤ tℓ =
    		 8C_,,∫_𝕊^2_+
    		U(v,t)
    		 S(v),

which we see is identical to (<ref>), and we have thus proved that X_M,R converges in distribution to X as we take M→∞ and then R→∞. This concludes the proof of Theorem <ref>. 

Computing the mean value 

We will determine the mean value (<ref>) of X; to do this we
exploit the integral geometry interpretation of the random
variable X.  By symmetry it suffices to restrict to directed lines
ℓ with v∈𝕊^2_+.  For fixed v∈𝕊^2_+,
denote by Q(v)=(K+span(v))∩ v^ the set of
q∈ v^ such that the directed line ℓ parametrized by
(v,q) intersects K.  We note that X A(q) is a volume element
of the box K for any fixed v∈𝕊^2_+, and thus
integrating X A(q) over all q yields the volume of the box.
Hence the mean value is 

    X=
    	8C_,,∫_𝕊^2_+∫_Q(v)X A(q) S(v) =
    	8C_,,∫_𝕊^2_+ S(v) =
    	2++.

In n dimensions we get a normalizing factor
Area(K)/2· 2^n∫_𝕊^n-1_+ v_n
S(v), so with the aid of the Lemma <ref> in the
Appendix, it follows that the mean value in n dimensions is 

    X =
      12^n   1π|𝕊^n|2^nArea(K)/2
      2^nVol(K)|𝕊^n-1|2^n = 2π|𝕊^n-1||𝕊^n|Vol(K)Area(K)

where Area(K) is the (n-1)-dimensional surface area
of the box K, and Vol(K) is the volume of the box
K. 



§ PROOF OF THEOREM <REF>

Using formula (<ref>) in dimension n=2, we
get

    cdf_X(t) =
        1-
        ∫_
          v∈𝕊^1_+
    ≤/t 
    ≤/t(
               (-t) +
               (-t)
            )
           S(v)
        ∫_𝕊^1_+( /+/)
           S(v) .


We use polar coordinates =cosθ, =sinθ so that S(v)=θ. Then the above becomes 

    1-
        ∫_cos (min(a/t, 1))  ^ sin (min(b/t, 1)) (
              cosθ+sinθ-2tsinθcosθ)
          θ∫_0^π/2( cosθ + sinθ)
        θ = 
    
        1-
        1/+[
            sinθ-cosθ+tcos^2θ]
          _cos (min(a/t, 1))  ^ sin (min(b/t, 1)) .

    

The numerator of the second term may be written

    χ(< t)( ·/t-√(1-^2/t^2)+t(1-^2/t^2) ) + 
      χ(≥ t)(  - · 0 + t· 0  ) + 
    
      -χ(< t)( √(1-^2/t^2) - ·/t + t·^2/t^2)
      -χ(≥ t)( · 0 - +t )

which can be simplified to
 
    χ(< t)( t--√(1-^2/t^2)) + 
      χ(< t)( t - -√(1-^2/t^2)) + 
      (  +-t ).

Inserting this into (<ref>) and differentiating yields Theorem <ref>.



§ PROOF OF THEOREM <REF>

We will evaluate the cumulative distribution function (<ref>) and then differentiate.
The denominator of the second term of (<ref>) is

    ∫_𝕊^2_+(++) S(v) =
            π4 (++),

as may be evaluated by switching to spherical coordinates.
Define

    f(,,)  S(v),
    
    	g(,,)  -2t S(v),
    
    	h(,,)  3t^2  S(v).

By symmetry, we have

    f(,,) =
    	 S(v) =
    	 S(v),
    
    	f(,,) =
    	 S(v) =
    	 S(v),
    
    	g(,,) =
    	-2t S(v) =
    	-2t S(v), 
    
    	g(,,) =
    	-2t S(v) =
    	-2t S(v),

and thus we can write the numerator in the second term of (<ref>) as

    f(,,)+
    	f(,,)+
    	f(,,)+
    	g(,,)+
    	g(,,)+
    	g(,,)+
    	h(,,).

Exploiting the symmetries, it suffices to evaluate h(,,), g(,,) and f(,,) (note the order of the arguments to f). 
We will evaluate these integrals by switching to spherical coordinates, but first we need to parametrize the part of the sphere inside the box 0≤≤/t, 0≤≤/t, 0≤≤/t.






 
  Fix t∈(0,√(^2+^2+^2)). We have
  
    F(,,)  S(v) = 
    (
            ∫_^θ_∫_0^π/2 +
            ∫_θ_^∫_φ_^π/2 -
            ∫_θ_^∫_φ_^π/2) F̃(θ,φ) sinθφθ

  for any integrable function F: 𝕊^2_+→, where F̃(θ,φ)  F(sinθcosφ,
  sinθsinφ, cosθ), where 
  
    /t, 
    θ_   max(, /t), 
    θ_   max(, /t), 
       √(^2+^2)t, 
    φ_   /tsinθ    (whenever ≤ tsinθ),
    φ_   /tsinθ    (whenever ≤ tsinθ).

  and where we have used the shorthand umin(u,1).



We will parametrize the set of points v=(,,) on the sphere
𝕊^2 such that 

    0< ≤/t,
    
    	0< ≤/t, 
    
    	0< ≤/t.

Switch to spherical coordinates =sinθcosφ, =sinθsinφ, =cosθ. The non-negativity conditions of (<ref>) are equivalent to the condition θ,φ∈(0,π/2). For such angles, the condition ≤/t is equivalent to

    /t≤θ,

and the conditions ≤/t, ≤/t are equivalent to

    tsinθ≤φ≤tsinθ.

The interval (<ref>) is non-empty for precisely those θ∈(0,π/2) such that θ≤ since

    1 ≤tsinθ^2 + tsinθ^2 
      1 ≤(tsinθ)^2 + (tsinθ)^2 
    sinθ≤√(^2+^2)t  θ≤√(^2+^2)t.

Thus we may restrict θ to the interval given by the inequalities

    ≤θ≤.

Note that we have ≤ for all t≤√(^2+^2+^2) since

    ≤
      1 ≤/t^2 + √(^2+^2)/t^2 
    
      1 ≤(/t)^2 + (√(^2+^2)/t)^2 
      t^2 ≤^2+^2+^2.

We conclude that we can write

    F(,,)  S(v) =
    	∫_^∫_/tsinθ^/tsinθF̃(θ,φ) sinθφθ.


For θ∈(0,π/2), note that /tsinθ is defined  precisely when
/t≤θ
and that /tsinθ is defined precisely when
/t≤θ.
  We have < θ_ if and only if t< √(^2+^2), and we have <θ_ if and only if t<√(^2+^2). Moreover we note that we always have θ_, θ_∈ [, ].

Let us rewrite the integration limits in the right-hand side of (<ref>) in terms of φ_ and φ_.
A priori, we need to distinguish between the two cases θ_≤θ_ and θ_< θ_. If θ_≤θ_ then we get


    (∫_^∫_x/tsinθ^y/tsinθ) =
    	(
    		∫_^θ_∫_0^π/2 +
    		∫_θ_^θ_∫_φ_^π/2 +
    		∫_θ_^∫_φ_^φ_) = 
    (
    		∫_^θ_∫_0^π/2 +
    		∫_θ_^∫_φ_^π/2 -
    		∫_θ_^∫_φ_^π/2 +
    		∫_θ_^∫_φ_^π/2 -
    		∫_θ_^∫_φ_^π/2) = 
    (
    		∫_^θ_∫_0^π/2 +
    		∫_θ_^∫_φ_^π/2 -
    		∫_θ_^∫_φ_^π/2).

If on the other hand θ_<θ_ then 

    (∫_^∫_x/tsinθ^y/tsinθ) =
    	(
    		∫_^θ_∫_0^π/2 +
    		∫_θ_^θ_∫_0^φ_ +
    		∫_θ_^∫_φ_^φ_) = 
    (
    		∫_^θ_∫_0^π/2 +
    		∫_θ_^θ_∫_0^π/2 -
    		∫_θ_^θ_∫_φ_^π/2 +
    		∫_θ_^∫_φ_^π/2 -
    		∫_θ_^∫_φ_^π/2)

which we see is identical to (<ref>). Combining (<ref>) and (<ref>) we get the conclusion of the lemma.


Applying Lemma <ref> 
we get

    h(,,) = 3t^2  S(v) = 
    
    	3t^2
    	(
    		∫_^θ_∫_0^π/2 +
    		∫_θ_^∫_φ_^π/2 -
    		∫_θ_^∫_φ_^π/2)
    		(sin^2θcosθcosφsinφ) sinθφθ.

An antiderivative of the integrand cosφsinφ·sin^3θcosθ with respect to φ is
-1/2cos ^2φsin ^3θcosθ, and thus the above is

    3t^2
    	(
    		∫_^θ_. cos ^2φ|_φ=0 +
    		∫_θ_^. cos ^2φ|_φ=φ_ -
    		∫_θ_^. cos ^2φ|_φ=φ_)
    		1/2sin ^3θcosθθ = 
    
    	3t^2
    	(
    		∫_^θ_
    			1 +
    		∫_θ_^a^2t^2sin^2θ +
    		∫_θ_^(b^2t^2sin^2θ-1)
    	)
    		1/2sin ^3θcosθθ = 
    3/2(
    		∫_^θ_
    			t^2sin ^3θcosθθ +
    		∫_θ_^
    			a^2 sinθcosθθ +
    		∫_θ_^(b^2sinθ - t^2sin^3θ) cosθθ) = 
    3/2(
    		[
    			t^21/4sin ^4θ]_^θ_ +
    		[
    			a^2 1/2sin^2 θ]_θ_^+
    		[
    			b^21/2sin^2θ-t^21/4sin^4θ]_θ_^).


Next consider

    g(,,) = -2t S(v) = 
    
      -2t(
        ∫_^θ_∫_0^π/2 +
        ∫_θ_^∫_φ_^π/2 -
        ∫_θ_^∫_φ_^π/2)
        (sin^2θcosφsinφ) sinθφθ.

An antiderivative of the integrand cosφsinφ·sin^3θ with respect to φ is
-1/2cos ^2φsin ^3θ, and thus the above is

    g(,,) = -2t S(v) = 
    
      -t(
        ∫_^θ_.cos ^2φ|_φ=0 +
        ∫_θ_^.cos ^2φ|_φ=φ_ -
        ∫_θ_^.cos ^2φ|_φ= φ_)
        sin ^3θθ = 
    
      -t(
        ∫_^θ_ 1 +
        ∫_θ_^a^2t^2sin^2θ +
        ∫_θ_^(b^2t^2sin^2θ-1)
      )
        sin ^3θθ = 
    
      -t(
        ∫_^θ_sin ^3θθ +
        ∫_θ_^a^2sinθt^2θ +
        ∫_θ_^(b^2sinθt^2-sin ^3θ)θ) = 
    
      -t(
        [cos^3θ/3-cosθ]_^θ_  +
        a^2t^2[-cosθ]_θ_^  +
        [-b^2cosθt^2-cos^3θ/3+cosθ]_θ_^).

We obtain g(,,) and g(,,) by switching the roles of
,, in (<ref>). We remark that trying to obtain
g(,,) and g(,,) directly, by integrating  and
, respectively, by first integrating with respect to
φ, taking the limits φ→φ_a and
φ→φ_b, and then finding an antiderivative with respect
to θ, seem to result in much more complicated expressions.

Finally consider

    f(,,) =
    	 S(v) = 
    (
    		∫_^θ_∫_0^π/2 +
    		∫_θ_^∫_φ_^π/2 -
    		∫_θ_^∫_φ_^π/2)
    		(sinθsinφ) sinθφθ.

An antiderivative of the integrand sinφ·sin^2θ with respect to φ is -cosφ·sin^2θ, and thus the above is

    (
    		∫_^θ_. cosφ|_φ=0 +
    		∫_θ_^. cosφ|_φ=φ_ -
    		∫_θ_^. cosφ|_φ=φ_)
    		sin ^2θθ = 
    (
    		∫_^θ_
    			1 +
    		∫_θ_^atsinθ -
    		∫_θ_^√(1-b^2t^2sin^2θ))
    		sin ^2θθ = 
    (
    		∫_^θ_sin^2θθ +
    		∫_θ_^asinθtθ -
    		∫_θ_^√(sin^2θ-b^2t^2)sinθθ)
    	= 
    (
    		1/2[
    			θ-sinθcosθ]_^θ_ +
    		[
    			-acosθt]_θ_^ -
    		∫_θ_^√(1-b^2t^2 - cos^2θ)sinθθ)

where the last integral inside the parentheses may be written as

    [
    			-1/2(cosθ√(1-b^2t^2-cos^2θ)+(1-b^2t^2) (cosθ/√(1-b^2t^2-cos^2θ)))
    		]_θ_^ = 
    [
    			-1/2(cosθ√(sin^2θ-b^2t^2)+(1-b^2t^2) (cosθ/√(sin^2θ-b^2t^2)))
    		]_θ_^

whenever θ_<π/2, by using the fact that 1/2(x √(c-x^2)+c (x/√(c-x^2))) is an antiderivative of √(c-x^2) with respect to x when c is a constant.

We obtain f(,,) and f(,,) by switching the roles of ,, in (<ref>).

It remains to insert the limits ,θ_,θ_, into the antiderivatives (<ref>), (<ref>) and (<ref>) above. Noting that ,θ_,θ_, are expressed in terms of piecewise-defined functions, the following manipulations will be useful.
For any function ψ, we have

    ψ()
      =ψ(cos/t)χ_ + ψ(cos 1)(1-χ_) 
    
      =(ψ(cos/t)-ψ(0))χ_ + ψ(0)

where χ_χ(t>). Similarly,

    ψ()
      =(ψ(sin√(^2+^2)/t)-ψ(π/2))χ_, +
       ψ(π/2)

where χ_,χ(√(^2+^2)>t), and

    ψ(θ_)=
      (1-χ_)ψ(π/2) +
      (χ_-χ_,) ψ(sin/t) +
      χ_,ψ(cos/t) 
    
      = χ_,·(
        ψ(cos/t)-ψ(sin/t)
      )+
      χ_·(
        ψ(sin/t)-ψ(π/2)
      )+
      ψ(π/2)

and similarly, ψ(θ_) can be written as

    χ_,·(
        ψ(cos/t)-ψ(sin/t)
      )+
      χ_·(
        ψ(sin/t)-ψ(π/2)
      )+
      ψ(π/2).


With this we can evaluate [ψ]_^θ_, [ψ]_θ_^, [ψ]_θ_^. But since we know that we will get a function symmetric with respect to the values ,,, it suffices to keep only those terms with χ_ and χ_,, say, and then the other terms may be evaluated by just switching the order of ,,. 
Upon inserting the limits and differentiating, one obtains (after tedious calculations) that 
  
    pdf_X(t) = F(,,,t) + F(,,,t) + F(,,,t)3 π  t^3 (++)

  where
  
    F(,,,t)     (8 t^3-3t^4) + 
        χ(t≥)((6t^4-^4+6 π^2 )-   (8  t^3-3t^4)-4(+)√(t^2-^2)(^2+2t^2))+
    
    
        χ(t≥√(^2+^2))[
                  ^4+^4-9 t^4-6 ^2 ^2+√(t^2-^2-^2) 4( ^2 + ^2 +2 t^2)+
            4√(t^2-^2)(^2+2t^2)-12 ^2 ·arctan(√(t^2-^2-^2))+
            4√(t^2-^2)(^2+2t^2)-12 ^2 ·arctan(√( t^2-^2-^2)/) ].

Rewriting F as a piecewise function, we get Theorem (<ref>).



§ PROOF OF THEOREM <REF>

Consider the distribution of the random variable
Y_M,N.  Since we record the same number of bounces
for each choice of angle φ we may replace the M-particle system
with a one particle system Y_N as follows: randomly select, with
uniform distribution, the angle φ and generate N bounce lengths
and randomly select one of these bounce lengths (with uniform
distribution); by the strong law of large numbers, Y_M,N converges
in distribution 
to Y_N as M→∞.

We now determine the limit distribution of Y_N.

As before, we first unfold the motion, and replace motion in a box
with specular reflections on the walls with motion in ^2; see
Figure <ref>.  The path lengths between bounces is then
the same as the lengths between the intersections with horizontal or
vertical grid lines.  To understand the spatial distribution, we
project the dynamics to the torus ^2/Λ where Λ is
the lattice

    Λ = {  (n_1,n_2) : n_1,n_2 ∈},

and we may identify the torus with the rectangle [0,]×[0,].


Let us first consider the motion of a single particle with an arbitrary
initial position, and direction of motion given by an angle φ.
Taking symmetries into account, we may assume
that φ∈ [0,π/2].  (Note that d φ/π/2
gives a probability measure on these angles.)

If the particle travels a large distance R>0, the number of intersections with
horizontal, respectively vertical, grid lines is R sinφ/ +O(1), respectively R cosφ/ +O(1).

Thus, in the limit R →∞, the probability of a line segment
beginning at a horizontal (respectively vertical) grid line is given
by P_h, respectively P_v (here we suppress the dependence on
φ) where

    P_hsinφ//sinφ/+cosφ/,   
    P_vcosφ//sinφ/+cosφ/.

The unfolded flow on the torus is ergodic for almost all φ,
and thus the starting points of the line segments becomes uniformly
distributed as R→∞ for almost all φ.

Let 

    T = T(φ) /cosφ.

Since sinφ = √(T^2-^2)/T, we obtain that

    P_h = √(T^2-^2)/+√(T^2-^2), 
      
    P_v = /+√(T^2-^2).



Let θ = arctan/ denote the angle of the diagonal in the
box, and assume that 0 ≤φ≤θ.  We then
observe the following regarding the line segment lengths.

First, if the segment begins at a horizontal line, it must end at a
vertical line, and the possible lengths of these segment lie between
0 and T.  We find that these lengths are uniformly distributed in
[0,T] since the starting points of the segments are uniformly
distributed.

On the other hand, if the line segment begins at a vertical line, it
can either end at a vertical or horizontal line.  Since the
starting points are uniformly distributed, the former happens with probability

    tanφ/ = 
    √(T^2-^2)//  = 
    √(T^2-^2)/

and the length of the segment is again uniformly distributed in
[0,T], whereas the latter happens with probability 

    - tanφ/ = 
    1- √(T^2-^2)/

in which case the segment is always of length T.  

Now, φ∈ [0,θ] implies that T ∈
[,√(^2+^2)], and noting that 

    d φ/ T = /T √(T^2-^2)

we find that the probability of observing a line segment of length t
is the sum of a “singular part” (the segment begins and ends on
vertical lines; note that all such segments have the same
lengths) and a “smooth part” (the segment does not begin and end on
vertical lines).  Moreover, the smooth part contribution equals

    1/π/2∫_max(,t)^√(^2+^2)1/T(
    P_h + P_vtanφ/) d φ/ T T

which, on inserting (<ref>), equals

    1/π/2∫_max(,t)^√(^2+^2)1/T·(
      √(T^2-^2)/+√(T^2-^2) + 
      /+√(T^2-^2)tanφ/) ·/T √(T^2-^2) T
      =
    1/π/2∫_max(,t)^√(^2+^2)1/T·(
      √(T^2-^2)/+√(T^2-^2) + 
      /+√(T^2-^2)√(T^2-^2)/) ·/T √(T^2-^2) T
      =
    1/π/2∫_max(,t)^√(^2+^2)2/+√(T^2-^2)· T/T^2 .



On the other hand, the “singular part contribution”, provided t
≥, to the probability of a segment having length t equals

    P_v/π/2·- tanφ/·d φ/dt 
    =
    1/π/2·/+√(t^2-^2)·(1- √(t^2-^2)/) ·/t √(t^2-^2)
    =
    1/π/2·/t  (+√(t^2-^2)) √(t^2-^2)·( - √(t^2-^2)).




In case θ≤φ≤π/2, a similar argument (we simple
reverse the roles of  and ) shows that the smooth contribution
equals

    1/π/2∫_max(,t)^√(^2+^2)2/+√(T^2-^2)· T/T^2

and that the singular contribution (if t ≥) equals

    1/π/2·/t  (+√(t^2-^2)) √(t^2-^2)·(- √(t^2-^2)) .




Thus, if we let P_sing(t) denote the “singular
contribution” to the probability density function we find the following: if t <, then

    P_sing(t) = 0

if t ∈ [,], then

    P_sing(t) = 
    1/π/2·(- √(t^2-^2)) /t  (+√(t^2-^2)) √(t^2-^2)

and if t ∈ [, √(^2+^2)], then

    P_sing(t)= 
    1/π/2·(
    (- √(t^2-^2))/t  (+√(t^2-^2)) √(t^2-^2) 
    +
    (- √(t^2-^2))/t  (+√(t^2-^2)) √(t^2-^2)).



  Note that P_sing has a singularity of type
  (t-)^-1/2 just to the right of t= (and similarly just to the right of t=).
  In a sense this singularity arises from
  the singularity in the change of variables φ↦ T since
  d φ/dT = /T √(T^2-^2).

  The reason for the singularities in the spreading model for n=2 is
  similar, as the spreading model can be obtained from the absorption model by a smooth
  change of the angular measure.



Similarly, the “smooth part” of the contribution is (for
t ∈ [0,√(^2+^2)]) given by

    P_smooth(t) =
    1/π/2(
    ∫_max(,t)^√(^2+^2)2/+√(T^2-^2)· T/T^2 
    +
    ∫_max(,t)^√(^2+^2)2/+√(T^2-^2)· T/T^2)





Hence the probability density function of the distribution of the segment length t is given by 

    pdf_Y(t) = P_sing(t) +  P_smooth(t).



We will now evaluate P_smooth(t).
An antiderivative of 2/+√(T^2-^2)·1/T^2 with respect to T for T∈(,√(^2+^2)) is

    2 ( √(T^2-^2)- )/T (^2+^2)+2 (tanh
       ^-1(T/√(^2+^2))-tanh ^-1(√(T^2-^2)√(^2+^2)/T ))/(^2+^2)^3/2

where tanh(z)=1/2log1+z/1-z for z<1. (A quick calculation shows that √(T^2-^2)√(^2+^2)/T  < 1 whenever <T<√(^2+^2).) We can rewrite (<ref>) as

    2 ( √(T^2-^2)- )/T (^2+^2)+log((√(^2+^2)+T) (T -√(T^2-^2)√(^2+^2))/(√(^2+^2)-T) (T  + √(T^2-^2)√(^2+^2)))
      /(^2+^2)^3/2

By l'Hôpital's rule we have

    lim_T→√(^2+^2)+T -√(T^2-^2)√(^2+^2)/√(^2+^2)-T  = 
      lim_T→√(^2+^2)+-T/√(T^2-^2)√(^2+^2)/-1  =
      
      ^2/

so the limit of (<ref>) as T→√(^2+^2)+ is

    log(
          (^2/)·(√(^2+^2)+√(^2+^2)) /(√(^2+^2)+√(^2+^2)))
      /(^2+^2)^3/2 =
      2log( /)
      /(^2+^2)^3/2.

The limit of (<ref>) as T→+ is

    -2 /(^2+^2)+2 tanh^-1(/√(^2+^2))/(^2+^2)^3/2.

Thus, assuming <, we can write π/2P_smooth(t) as

    2 (+)/(^2+^2)-2 /(^2+^2)^3/2(tanh^-1(/√(^2+^2)) + tanh^-1(/√(^2+^2)))

if t<,, or as

    2  +2 t - 2√(t^2-^2)/t (^2+^2)+
    2 (-tanh
       ^-1(t/√(^2+^2))+tanh ^-1(√(t^2-^2)√(^2+^2)/t ) - tanh^-1(/√(^2+^2)))/(^2+^2)^3/2

if <t< or as

    22- √(t^2-^2) - √(t^2-^2)/t (^2+^2)+ 
    2 (-2tanh
       ^-1(t/√(^2+^2))+tanh ^-1(√(t^2-^2)√(^2+^2)/t ) +
        tanh ^-1(√(t^2-^2)√(^2+^2)/t )
        )/(^2+^2)^3/2

if ,<t. Adding P_sing(t) to this, we get Theorem <ref>.








§ CALCULATION OF AN INTEGRAL



 
Write |𝕊^n-1| for the (n-1)-dimensional surface area of the sphere
𝕊^n-1⊆^n. Then we have 
  
    ∫_𝕊^n-1_+ v_n  S(v) = 1π|𝕊^n|2^n.

where 𝕊^n-1_+ 𝕊^n-1∩ (0,∞)^n is the part of the sphere 𝕊^n-1 with positive coordinates.


  We may parametrize v=(v_1,…,v_n)∈𝕊^n-1_+ with
  
    v_1    = cosθ_1 
    
        v_2    = sinθ_1 cosθ_2 
    
        v_3    = sinθ_1 sinθ_2 cosθ_3 
       ⋮
    
        v_n-1   = sinθ_1⋯sinθ_n-2cosθ_n-1
    
        v_n    = sinθ_1⋯sinθ_n-2sinθ_n-1

  for θ_1,…,θ_n-1∈(0,π/2). We have the spherical area element
  
    S(v) = sin^n-2θ_1 sin^n-3θ_2 ⋯sinθ_n-2θ_1⋯θ_n-1.

  Thus we get
  
    ∫_𝕊^n-1_+ v_n  S(v) = 
        ∏_i=1^n-1∫_0^π/2sin^n-1-iθ_i θ_i.

  Introducing an additional integration variable θ_n, we recognize the integrand as the spherical area element in n+1 dimensions, and thus the above is
  
    1∫_0^π/2θ_n∏_i=1^n∫_0^π/2sin^n-1-iθ_i θ_i =
        1π/2|𝕊^n|2^n+1.

  since ∫_𝕊^n_+ S(v) = |𝕊^n|/2^n+1.








99

barra-gaspard-level-spacings-quantum-graphs
F. Barra and P. Gaspard.
On the level spacing distribution in quantum graphs.
J. Statist. Phys., 101(1-2):283–319, 2000.



  
boldrighini-bunimovich-sinai-boltzmann-lorenz
C. Boldrighini, L. A. Bunimovich, and Y. G. Sinaĭ.
On the Boltzmann equation for the Lorentz gas.
J. Statist. Phys., 32(3):477–501, 1983.

bourgain-golse-wennberg-lorentz-gas
J. Bourgain, F. Golse, and B. Wennberg.
On the distribution of free path lengths for the periodic Lorentz
  gas.
Comm. Math. Phys., 190(3):491–508, 1998.

bunimovich-sinai
L. A. Bunimovich and Y. G. Sinaĭ.
Statistical properties of Lorentz gas with periodic configuration
  of scatterers.
Comm. Math. Phys., 78(4):479–497, 1980/81.


applied-paper M. Bäckström, S. Holmin, P. Kurlberg, D. Månsson.  Randomized Ray Tracing for Modeling UWB Transients in a Reverberation Chamber. In preparation.


chernov-mean-free-path-billiards
N. Chernov.
Entropy, Lyapunov exponents, and mean free path for billiards.
J. Statist. Phys., 88(1-2):1–29, 1997.

  
golse-wennberg-lorentz-gas
F. Golse and B. Wennberg.
On the distribution of free path lengths for the periodic Lorentz
  gas. II.
M2AN Math. Model. Numer. Anal., 34(6):1151–1163, 2000.

klain-rota-integral-geometry-book
D. A. Klain and G.-C. Rota.
Introduction to geometric probability.
Lezioni Lincee. [Lincei Lectures]. Cambridge University Press,
  Cambridge, 1997.


mansson-pers-comm
D. Månsson, personal communication.



marklof-strombergsson-annals-2010
J. Marklof and A. Strömbergsson.
The distribution of free path lengths in the periodic Lorentz gas
  and related lattice point problems.
Ann. of Math. (2), 172(3):1949–2033, 2010.

marklof-strombergsson-boltzman-grad-annals2011
J. Marklof and A. Strömbergsson.
The Boltzmann-Grad limit of the periodic Lorentz gas.
Ann. of Math. (2), 174(1):225–298, 2011.

marklof-strombergsson-gaps-between-logs
J. Marklof and A. Strömbergsson.
Gaps between logs.
Bull. Lond. Math. Soc., 45(6):1267–1280, 2013.

marklof-strombergsson-CMP-quasicrystal
J. Marklof and A. Strömbergsson.
Free path lengths in quasicrystals.
Comm. Math. Phys., 330(2):723–755, 2014.



santalo-book
L. A. Santaló.
Integral geometry and geometric probability.
Cambridge Mathematical Library. Cambridge University Press,
  Cambridge, second edition, 2004.
With a foreword by Mark Kac.

spohn-lorentz
H. Spohn.
The Lorentz process converges to a random flight process.
Comm. Math. Phys., 60(3):277–290, 1978.

wennberg-quasicrystal
B. Wennberg.
Free path lengths in quasi crystals.
J. Stat. Phys., 147(5):981–990, 2012.







