Linear Boltzmann dynamics in a strip with large reflective obstacles: stationary state and residence time

The presence of obstacles modify the way in which particles diffuse. In cells, for instance, it is observed that, due to the presence of macromolecules playing the role of obstacles, the mean square displacement ofbiomolecules scales as a power law with exponent smaller than one. On the other hand, different situations in grain and pedestrian dynamics in which the presence of an obstacle accelerate the dynamics are known. We focus on the time, called residence time, needed by particles to cross a strip assuming that the dynamics inside the strip follows the linear Boltzmann dynamics. We find that the residence time is not monotonic with the sizeand the location of the obstacles, since the obstacle can force those particles that eventually cross the strip to spend a smaller time in the strip itself. We focus on the case of a rectangular strip with two open sides and two reflective sides and we consider reflective obstaclea into the strip.


(Communicated by Laurent Desvillettes)
Abstract. The presence of obstacles modifies the way in which particles diffuse. In cells, for instance, it is observed that, due to the presence of macromolecules playing the role of obstacles, the mean-square displacement of biomolecules scales as a power law with exponent smaller than one. On the other hand, different situations in grain and pedestrian dynamics in which the presence of an obstacle accelerates the dynamics are known. We focus on the time, called the residence time, needed by particles to cross a strip assuming that the dynamics inside the strip follows the linear Boltzmann dynamics. We find that the residence time is not monotonic with respect to the size and the location of the obstacles, since the obstacle can force those particles that eventually cross the strip to spend a smaller time in the strip itself. We focus on the case of a rectangular strip with two open sides and two reflective sides and we consider reflective obstacles into the strip. We prove that the stationary state of the linear Boltzmann dynamics, in the diffusive regime, converges to the solution of the Laplace equation with Dirichlet boundary conditions on the open sides and homogeneous Neumann boundary conditions on the other sides and on the obstacle boundaries.
1. Introduction. Particle based studies of agent behavior can reproduce realistic scenarios. Many different real situations, such as grains discharging from a silo, traffic flow, and pedestrian dynamics, can be studied via particle based modeling. The focus, in this paper, is on the effect of obstacles on particle flows [15].
It is well known that the mean square distance traveled by particles undergoing Brownian motion is proportional to time. However, many experimental measures of molecular diffusion in cells show a sublinear behavior. This phenomenon, called anomalous diffusion, is in some cases explained as an effect of the presence of macromolecules playing the role of obstacles for diffusing smaller molecules [16,26,31,32] While in the case of anomalous diffusion the obstacles induce a sort of slowing down of the dynamics, on the contrary in many other different contexts it has been noticed and exploited the fact that the presence of an obstacle can surprisingly accelerate the dynamics improving, in particular, the flux rate of particles passing through a bottleneck.
In granular flows, for instance when grains discharge from a silo, it happens that the out-coming flow can be dramatically reduced due to clogging at the exit. In [33] it was proposed that this phenomenon is caused by the formation of arches. In the case of spherical grains it has been demonstrated that the presence of clogging in a three-dimensional silo critically depends on the ratio between the outlet size and the diameter of the particles [34]. A solution that is implemented to improve the granular flow is to place an obstacle above the silo exit [2,35] which prevents arches to be formed or to become stable.
A similar phenomenon is observed in pedestrian flows (see, e.g., [4,21,23,25] for reviews of models and related problems). In the case of pedestrian exiting a room under panic clogging at the door can be reduced by means of suitably positioned obstacles, see [1,Section 6.3] and [22,24]. It has also been noticed that a correct positioning can reduce injuries under panicked escape from a room thanks to the so called "waiting-room" effect [17]: pedestrians slow down and accumulate close to the obstacle so that the exit is decongested. We mention, here, that also the possibility of clustering far from the exit due to individual cooperation has been object of study in [10,13,14].
The a priori unexpected phenomena discussed above are a sort of inversion of the Braess' paradox [5,27] stating that adding a road link to a road network can cause cars to take longer to cross the network. Indeed, in the examples discussed above it seems that adding barriers results in a decrease of the time that particles need to cross a region of space. This is precisely the issue we address in this paper. Inspired by [11,12], we consider particles entering an horizontal strip through the left boundary and eventually exiting through the right one [20]. We compute the typical time needed to cross the strip, called the residence time, and analyze its dependence on the size and position of a reflecting obstacle positioned inside the strip [7,8]. Surprisingly, we find not monotonic behaviors of the residence time as a function of the side lengths of the obstacle and the coordinate of its center [9]. In particular, for suitably choices of the obstacles the residence time in presence of the barrier is smaller than the one measured for the empty strip. In other words, our results show that placing a suitable obstacle in the strip allows to select those particles that cross the strip in a shorter time. We also observe that the same obstacle, placed in different position, can either increase or decrease the residence time.
Inside the strip we consider particle moving according to the Markov process solving the linear Boltzmann equation. We stress that this is the first study of residence times by means of the methods of Kinetic Theory. In this case we calculate the residence time by directly simulating the motion of single particles by a Monte Carlo method. This dynamics should be consistent with the Lorentz process at appropriate regimes. The Lorentz gas model is a system of non-interacting particles moving in a region where static small disks are distributed according to a Poisson probability measure. This is a classical model for finite velocity random motions. Particles perform uniform linear motion up to the contact with a disk where they are elastically reflected. We study the system in a diffusive regime and we show that there exists a unique stationary solution which converges to the solution of a Laplace problem with mixed boundary conditions: Dirichlet boundary conditions on the vertical sides and homogeneous Neumann on the rest. We prove the convergence and check it numerically. Moreover, by constructing numerically the stationary profiles we qualitatively verify that the overall flux in presence of obstacles is decreasing, as expected by physical intuition. This holds even in those cases in which the residence time is smaller with respect to the empty strip case.
In Section 2 we introduce the model under investigation and we state our main theoretical results. In Section 3 we propose a Monte Carlo algorithm that we use to construct the stationary state of the linear Boltzmann model and we discuss the stationary profile in presence of large fixed obstacles. In Section 4 we discuss our results on residence time. Finally, Section 5 contains the proofs of the results we state in Section 2.
2. Model and results. We consider a system of light particles moving in the two-dimensional space. We choose as the domain a subset Ω of the finite strip (0, L 1 ) × (0, L 2 ) ⊂ R 2 . This strip has two open boundaries, that we think as the left side ∂Ω L = {0} × (0, L 2 ) and the right side ∂Ω R = {L 1 } × (0, L 2 ). The strip is in contact on the left side and on the right side with two mass reservoirs at equilibrium with particle mass densities ρ L and ρ R , respectively. Particles traveling into Ω are instead specularly reflected upon colliding with the upper side (0, L 1 ) × {L 2 } and lower side (0, L 1 ) × {0} of the strip.
We consider the case in which large fixed obstacles are placed in the strip so that the domain Ω is a connected set. These obstacles are convex sets with smooth reflective boundaries. We consider a generic configuration of a finite number of obstacles with positive mutual distance and positive distance from the sides of the strip. In the sequel we will call ∂Ω E the union of the obstacle boundaries and the upper and lower sides of the strip (see Figure 1). Therefore, when a particle reaches ∂Ω E it experiences a specular reflection.  The linear Boltzmann equation is a kinetic linear equation, combining free transport and scattering off of a medium. This equation consists of two terms: a free transport term and a collision operator L.
Let us consider the phase space Ω × S 1 , where S 1 := {v ∈ R 2 : |v| = 1}. We will consider the operator L with elastic collision kernel. The equation for (x, v) ∈ Ω×S 1 and positive times t reads where, by the elastic collision rule v = v − 2(n · v)n, the operator L is defined for any f ∈ L 1 (S 1 ) as Here n = n(δ) is the outward pointing normal to a circular scatterer of radius 1 at the point of collision among the particle with velocity v and the scatterer. So δ = sin α if α is the angle of incidence between v and n that has δ as impact parameter (see Figure 2); λ > 0 is a fixed parameter.
n v α  Figure 2. Elastic collision with a scatterers: impact parameter δ and angle of incidence α.
We denote by g ε the solution of the equation corresponding to the value of η ε , that is a positive parameter that we let go to +∞ as ε goes to 0 + . The choice of the kernel and the related parameters will be discussed at the end of this Section.
The equation describes the evolution of the density of particles, moving of linear motion and having random collisions, against a circular scatterer, that preserve the energy. The time between two consecutive jumps in the velocities is distributed with exponential law with mean value (λη ε ) −1 . Since both random collisions and hits against the elastic boundaries preserve the energy, the modulus of the velocity of a particle moving in Ω is constant, so we consider it to be equal to one.
On the elastic boundary ∂Ω E we impose reflective boundary condition and on the open boundary ∂Ω L ∪ ∂Ω R we set Dirichlet condition: where f B is defined in 4 below and v is given by the elastic collision rule v = v − 2(n · v)n. Here we denote by n = n(x) the inward pointing normal on the boundary ∂Ω of the domain. We consider as initial datum the function f 0 (x, v) ∈ L ∞ (Ω × S 1 ) and we define f B (not depending on t) as where 1/2π is the density of the uniform distribution on S 1 . We are interested in the study of the stationary problem associated to 1-3: We want to investigate the behavior of the solution g S ε of 5 and prove its convergence to the stationary solution of the diffusion problem in Ω with mixed boundary conditions given by Theorem 2.1. If ε > 0 is sufficiently small there exists a unique stationary solution as ε → 0, where ρ(x) is the solution to the problem 6. The convergence is in The choice of the elastic collision kernel for the operator L defined in 2 is due to the physical model we have in mind. We are considering a particle moving with initial velocity v ∈ S 1 and hitting an hard circular scatterer whose position is random. The random impact parameter δ chosen uniformly in [−1, 1] allows to individuate this collision. In a similar way we could let the particle move following the Lorentz process, that is moving freely in a region where static small disks of radius ε are distributed according to a Poisson probability measure and elastically colliding with those disks. In this case for suitable choices of the mean value of the Poisson distribution in terms of ε and η ε it could be possible to prove that the diffusive limit for the linear Boltzmann equation and the Lorentz process in a (small disks) low density limit are asymptotically equivalent in the limit ε → 0 (see [3] for the case of an infinite 2D slab with open boundary).
3. Numerical convergence of stationary linear Boltzmann. We investigate, here, the stationary solution of the linear Boltzmann equation from a numerical point of view. Our algorithm directly simulates the motion of single particles following the linear Boltzmann equation. In the simulations we exploit the interpretation of the linear Boltzmann equation as the equation describing a stochastic jump process in the velocities and we directly simulate the motion of single particles.
We will show that the numerical stationary solution that we construct is close to the solution of the associated Laplace problem 6 if the scale parameter ε is small enough, that is to say if the average time t m between two consecutive hits is sufficiently small. This time will be called in the sequel mean flight time.
We will construct the solution of the Laplace problem 6 in our geometry by using the COMSOL Multiphysics software.
We proceed in the following way: we consider particles entering in Ω from the reservoirs. A particle starts its trajectory from the left boundary ∂Ω L or from the right boundary ∂Ω R , where the mass density is ρ L and ρ R respectively. Therefore the number of particles we let enter from each side is chosen to be proportional to ρ L and ρ R . In other words, we select the starting side of the particle with probability ρ L /(ρ L + ρ R ) (left side) and ρ R /(ρ L + ρ R ) (right side) respectively. Then we draw uniformly the position x in ∂Ω L or ∂Ω R and the velocity v in S 1 with v · n(x) > 0, n(x) inward-pointing normal.
Once the particle started, it moves with uniform linear motion until it hits a scatterers or the elastic boundary ∂Ω E . We pick t h , the time until the hit with a scatterers, following the exponential law of mean t m , with t m a fixed parameter. The particle travels with velocity v for a time t h . If during this time it hits the elastic boundary ∂Ω E , its velocity changes performing an elastic collision. At time t h we simulate an hit with a scatterers by picking an impact parameter δ uniformly in [−1, 1] and changing the velocity from v to v = v − 2(v · n)n, where δ = sin α and n = n(δ) is the outward pointing normal to the scatterers such that the angle of incidence between n and v is α (see Figure 2).
We proceed as before by letting the particle move until it leaves Ω by reaching again the open boundary ∂Ω L ∪ ∂Ω R . Then the particle exits from the system and we are ready to simulate another particle. We simulate a number N of particles.
The random number generator we use in our simulation is the Mersenne Twister [29,30].
We want to construct the stationary solution of equation 1-3. Note that we can simulate particles one by one since in the considered model particles are not interacting. Moreover, being the stationary state not dependent on the initial datum f 0 , we consider in this algorithm only particles starting from the reservoirs. We assume that in the stationary state the density of particles in a region is proportional to the total time spent by all particles in that region. Moreover, due to isotropy, there is no preferential direction for the velocity, so the stationary g S ε is not dependent on v.
We divide the space Ω in equal small square cells. For every particle we calculate the time that it spends in every cell. Then we calculate the total time that particles spend in each cell. For t m going to zero and N very big, in every infinitesimal region of Ω the stationary density has to be proportional to the total time spent by particles in that region. Considering sufficiently small cells, for t m small enough and a number of particles simulated N big enough, the total time spent in the cell we calculate is proportional to our numerical stationary solution.
We construct with our algorithm a grid of sojourn times in the cells. The last step we have to do is to normalize it. It is sufficient to multiply by a constant, obtained by imposing the correct value of g S ε in a point (e.g., the boundary datum). We call our numerical solution h tm (x). So we fix a cell in contact with the reservoir where we calculated a sojourn time t c and consider the value f B of the stationary solution. We choose as multiplication constant c = f B (x)/t c . So the simulated solution h tm (x) is constructed by multiplying the sojourn time in each cell for this constant c.
In the sequel we will show that for t m sufficiently small and N big enough, the simulated solution h tm (x) well approximates the solution ρ(x) of the associated Laplace problem.
All the simulations we are going to discuss in this Section are performed with N = 5 · 10 7 particles.
Let us preliminary consider the case of Ω = (0, 4) × (0, 1) in absence of obstacles. We fix mass densities at the reservoirs ρ L = 1 and ρ R = 0.5. In this first case there is no dependence on the vertical coordinate in the solution of the Laplace problem 6. Indeed we know that the problem has analytic solution ρ(x 1 , x 2 ) = 1 − x 1 /8, where we are denoting by (x 1 , x 2 ) ∈ R × R the spatial coordinates x in Ω. We divide the domain in 200 × 50 equal square cells and we consider simulations with different values of t m , to understand which values of the mean time t m provide a good approximation of the solution we are looking for.
In Figures 3 and 4 we show that the choice of t m of the order of 10 −2 is suitable for our purpose. Indeed in Figures 3 we compare the simulated solution h tm with the analytic solution for the values t m = 2 · 10 −1 , t m = 10 −1 , t m = 2 · 10 −2 . We see in a 3D plot and in a 2D plot, obtained from the previous one by averaging on the x 2 variable, that h tm becomes closer to ρ(x 1 ) when t m decreases. So in Figures 4 we fix the parameter t m = 10 −2 and we verify that h tm is close to the function ρ(x 1 ) showing the relative error |h tm − ρ|/ρ.  We consider now the more interesting case with presence of obstacles in the strip. Our domain Ω is the strip (0, 4) × (0, 1) minus the obstacles. We fix again mass densities at the reservoirs ρ L = 1 and ρ R = 0.5. We fix again t m = 10 −2 , since we have shown that in the empty case this choice for the exponential clock allows to construct a numerical solution h tm that is close to the analytical solution ρ(x 1 , x 2 ) of the associated Laplace problem 6. We propose different situations for the domain Ω and we show in Figures 5 -7 that in each case the simulation algorithm works correctly. We compare our numerical solution with the solution ρ(x 1 , x 2 ) of the associated Laplace problem 6. We show the plots of the h tm and ρ and the map of the relative error |h tm − ρ|/ρ as in Figure 4 .
The first case we consider is the presence of a big squared obstacle with side 8 · 10 −1 , in different positions into the strip. In Figure 5 the results on two different positions are presented.  Another interesting case is the presence of a very thin and tall obstacle placed vertically inside the strip. We show it in Figure 6, by picking a thin obstacle of height of 8 · 10 −1 .
The last case we want to present is the presence in the strip of two obstacles. In Figure 7 we consider two different situations: in the first we set in the strip two squared obstacles with sides 6 · 10 −1 long, in the second we place into the strip two rectangular obstacles of sides 4 · 10 −1 and 7 · 10 −1 .
Note that due to the presence of obstacles the solutions are not independent of the vertical coordinate x 2 anymore as it was in the empty strip case. However, we can notice that before and beyond the obstacles in the x 1 direction the stationary states are closer to a flat state than in the empty case, with a steeper slope in the tight channels at sides of the obstacles. The total stationary mass flux through any vertical line {x 1 } × (0, 1) ∩ Ω does not depend on x 1 . Indeed, this should follow from the Fick's law, that we expect to be valid also in presence of obstacles (in absence   of obstacle, being the limiting problem one-dimensional, the Fick's law holds as shown in [3]), together with the divergence theorem and the fact that the boundary conditions are homogeneous on ∂Ω E . The Fick's law would tell us that, in presence of obstacles, the total flux on the lines {x 1 = a} ∩ Ω is smaller than in the empty case, as it is possible to see focusing on the vertical lines before the obstacles. In this sense, and opposite to what happen in the case of the study of the residence time, we find on the flux the intuitive result we expected.

Residence time.
We consider the domain (0, L 1 )×(0, L 2 ) with the same boundary conditions as in Section 2, namely, reflecting horizontal boundaries and open vertical boundaries. As before Ω denotes a subset of this domain obtained by placing large fixed reflecting obstacles. Particles in Ω move according to the Markov process solving the linear Boltzmann equation and described in detail in Section 3.
In Section 3 we investigated the stationary state of the system and we demonstrated that, provided the mean flight time t m is sufficiently small, the stationary state is very well approximated by the solution of the Laplace problem 6 even in presence of obstacles. We have also noted that, due to the presence of obstacles, the total flux crossing the strip is smaller with respect to the one measured in absence of obstacles. This implies that if we consider a fixed number of particles entering the strip throught the left boundary, the number of them exiting through the right boundary decreses when an obstacle is inserted. In our simulations we remark that the ratio between the number of particles exiting through the left boundary in presence of an obstacle and in the empty strip case does not depend very much on the geometry of the obstacle and, in the worst case we considered, it is approximatively equal to 1/5. Detailed data for the different cases we studied are reported in the figure captions of this section.
In this section, on the other hand, we focus on those particles that do the entire trip, that is to say they enter through the left boundary and eventually exit the strip through the right one. Limiting our numerical computation to these particles, we measure the average time needed to cross the strip, also called the residence time and discuss its dependence on the size and on the position of a large fixed obstacle placed in the strip. The surprising result is that the residence time is not monotonic with respect to the obstacle parameters, such as position and size. More precisely, we show that obstacles can increase or decrease the residence time with respect to the empty strip case depending on their side lengths and on their position. Moreover, in some cases, by varying only one of these parameters a transition from the increasing effect to the decreasing effect is observed.
In some cases we observe that the residence time measured in presence of an obstacle is smaller than the one measured for the empty strip. In other words, we find that the obstacle is able to select those particles that cross the strip in a smaller time. More precisely, particles that succeed to cross the strip do it faster than they would in absence of obstacles.
We now discuss the different cases we considered. All numerical details are in the figure captions. The statistical error is not represented in the pictures since it is negligible and it could not be appreciated in the graphs. In each figure we draw a graph reporting the numerical data and a schematic picture illustrating the performed experiment. We first describe our result and at the end of this section we propose a possible interpretation.
In Figure 8 we report the residence time as a function of the obstacle height. The obstacle is placed at the center of the strip and its width is very small on the left and larger on the right. We notice that in the case of a thin barrier, the residence time increases with the height of the obstacle. On the other hand, for a wider obstacle, we do not find this a priori intuitive result, but we observe a not monotonic dependence of the residence time on the obstacle height. In particular, it is interesting to remark that if the obstacle height is chosen smaller that 0.65 the residence time is smaller than the one measured for the empty strip. This effect is even stronger if the width of the obstacle is increased (Figure 9).   Figure 8. Residence time vs. height of a centered rectangular obstacle with fixed width 4 · 10 −2 (on the left) and 4 · 10 −1 (on the right). Simulation parameters: L 1 = 4, L 2 = 1, t m = 2 · 10 −2 , total number of inserted particles 10 8 , the total number of particles exiting through the right boundary varies from 5.3 · 10 5 to 3.6 · 10 5 (on the left) and from 5.3 · 10 5 to 2.1 · 10 5 (on the right) depending on the obstacle height. The solid lines represent the value of the residence time measured for the empty strip (no obstacle).
In the left panel of Figure 10 we report the residence time as a function of the obstacle width. The obstacle is placed at the center of the strip and its height is fixed to 0.8. When the barrier is thin the residence time is larger than the one measured in the empty strip case. But, when the width is increased, the residence time decreases and at about 0.7 it becomes smaller than the empty case value. The minimum is reached at about 2.3, then the residence time starts to increase and when the width of the obstacles equals that of the strip the residence time becomes equal to the empty strip value. This last fact is rather obvious, indeed, in this case the strip reduces to two independent channels having the same width of the original strip.
In the right panel of Figure 10 a centered square obstacle is considered. We note that the residence time happens to be a monotonic decreasing function of the obstacle side length.
In Figure 11 we show that, and this is really surprising, the residence time is not monotonic even as a function of the position of the center of the obstacle. In the left panel a squared obstacle of side length 0.8 is considered, whereas in the right panel a thin rectangular obstacle 0.04 × 0.8 is placed in the strip. In both cases the residence time is not monotonic and attains its minimum value when the obstacle is placed in the center of the strip. It is worth noting, that in the case on the left when the position of the center lays between 1.5 and 2.5 the residence time in presence of the obstacles is smaller than the corresponding value for the empty strip.   Figure 9. Residence time vs. height of a centered rectangular obstacle with fixed width 8 · 10 −1 (on the left) and 12 · 10 −1 (on the right). Simulation parameters: L 1 = 4, L 2 = 1, t m = 2 · 10 −2 , total number of inserted particles 10 8 , the total number of particles exiting through the right boundary varies from 5.2 · 10 5 to 1.4 · 10 5 (on the left) and from 5.2 · 10 5 to 1.1 · 10 5 (on the right) depending on the obstacle height. The solid lines represent the value of the residence time measured for the empty strip (no obstacle).
Summarizing, the numerical experiments reported in Figures 8-11 show that the residence time strongly depends on the obstacle geometry and position. In particular it is seen that large centered obstacles favor the selection of particles crossing the strip faster than in the empty strip case.
A possible interpretation of these results can be given. The strip (0, L 1 ) × (0, L 2 ) is partitioned in the three rectangles L (the part on the left of the obstacles), R (the part on the right of the obstacles), and C = (0, L 1 ) × (0, L 2 ) \ (L ∪ R). The phenomenon we reported above can be explained as a consequence of two competing effects: the total time spent by a particle in the channels between the obstacle and the horizontal boundaries is smaller with respect to the time typically spent in C in the empty strip case because of the volume reduction due to the presence of the obstacle. On the other hand the times spent in L and in R are larger if compared to the times spent there by a particle in the empty strip case, due to the fact that it is more difficult to leave these regions and enter in the channels flanking the obstacle. The increase or the decrease of the residence time compared to the empty strip case depends on which of the two effects dominates the particle dynamics.
In Figure 12 we consider the geometry in the right panel of Figure 8. We compute the average time spent by particles in small squared cells (0.02 × 0.02). This local residence time in presence of the obstacle is larger than the one measured in the empty strip, indeed the gray surface in the picture is always above the black one. But, if the total residence time spent in the regions L, C, and R is computed, one discovers that the time spent in the region C decreases in presence of the obstacle,   whereas the time spent in L and R increases. Note that the local residence time in the cells belonging to the channels in C is larger with respect to the empty strip case, but the total time in C is smaller due to the fact that the available volume in C is decreased by the presence of the obstacle. Hence, the result in the right panel in Figure 8 can be explained as follows: if the height of the obstacle is smaller than 0.6 the effect in C dominates the one in L and R so that the total residence time decreases. On the other hand, when the height is larger than 0.6 the increase of the residence time in L and R dominates its decrease in C, so that the total residence time increases. The Figure 13, referring to the geometry in the left panel in Figure 10, and the Figure 14, referring to the geometry in the left panel in Figure 11, can be discussed similarly. We just note that in Figures 12 and 13 the circles and triangles, which correspond to the residence time in L and R, coincide due to the symmetry of the system. Indeed, in both cases the center of the obstacle is at the center of the strip.

5.
Proof of results. We prove Theorems 2.1 and 2.2. We firstly construct the solution of the linear Boltzmann problem in form of a Dyson series. Then we are able to prove the existence and uniqueness of the associated stationary problem. To do this we exploit the diffusive limit of the linear Boltzmann equation in a L ∞ setting and in a bigger domain containing Ω, by means of the Hilbert expansion method (see [3,6,18]). The stationary solution is constructed in the form of a      Figure 13. As in Figure 12 for the geometry in the left panel in Figure 10. In the left panel the width of the obstacle is 2.28.  Figure 12 for the geometry in the left panel in Figure 11. In the left panel the position of the center of the obstacle is 0.8.
Neumann series to avoid the exchange of the limits t → ∞, ε → 0, following the idea of [3]. Eventually we prove the convergence of the stationary state to the solution of the mixed Laplace problem. This also requires the Hilbert expansion method. The auxiliary results stated are proved after the main theorems. Let us consider the problem 1-3 with the datum g ε (x, v, 0) = f 0 (x, v) ∈ L ∞ (Ω × S 1 ). We can express the operator defined in 2 as (8) and I is the identity. Therefore the equation 1 can be written as We want to exploit the Duhamel's principle and express the solution as a series expansion. We consider the semigroup generated by A = (v · ∇ x + 2η ε λI). We recall that in the whole plane R 2 this semigroup acts as e −tA f (x, v) = e −2ληεt f (x − vt, v), while the semigroup generated only by the transport term v · ∇ x would be We want to consider the semigroup generated by A on our domain Ω initial datum f 0 and boundary conditions 3. Recall that ∂Ω E is a specular reflective boundary while on ∂Ω L ∪ ∂Ω R the system is in contact with reservoirs with particle densities f B (x, v). Since the equation describes the evolution of a particle moving in the space with velocity of modulus one, having random collisions with impact parameter δ, for any sequence of collision times and impact parameters t i , δ i we can construct the backward trajectory of a particle as long as it stays in Ω. Indeed the backward trajectory for a particle in (x, v) at time t starts by letting the particle move with velocity −v. For a time t − t 1 it does not hit any scatterers, but if the particle reaches the elastic boundary ∂Ω E during this time, the velocity −v becomes −v following the elastic collision rule −v = −(v − 2(n · v)n), where n is the inward pointing normal to Ω. After a time t−t 1 the particle performs a collision with impact parameter δ 1 that produces the velocity −v 1 . Then again the particles travels for a time t 1 − t 2 elastically colliding if touching the boundary ∂Ω E and so on until it reaches a reservoir or it has traveled for a total time t.
In the same way, given the sequence x, v, t 1 , . . ., t m , δ 1 , . . ., δ m , we define the flow Φ −t m (x, v, t 1 , . . . , t m , δ 1 , . . . , δ m ) as the backward trajectory starting from x with velocity v and having m transition in velocity obtained after a time t − t 1 , . . ., t i − t i+1 , . . ., t m (i = 1, . . . , m − 1) by a scattering with an hard disk with impact parameter respectively δ i (i = 1, . . . , m). We impose that the trajectories described by this flow Φ −t (x, v) make a change of velocity from v to v = v − 2(n · v)n any time the elastic boundary ∂Ω E is reached.
We define the function τ = τ (x, v, t, t 1 , . . . , t m , δ 1 , . . . , δ m ) that represents the time when the particle that is in (x, v) at time t leaves a reservoirs and it enters into the strip. So if the backward trajectory having collision times and parameters t 1 , . . . , t m , δ 1 , . . . , δ m reaches the boundary ∂Ω L ∪ ∂Ω R in the time interval [0, t], then it happens after a backward time t − τ . If the trajectory Φ −s (x, v) never hits the boundary ∂Ω L ∪ ∂Ω R for any time s ∈ [0, t] we set τ = 0.
We are now able to write the solution g ε (x, v, t) using the Duhamel's principle. The semigroup generated by A on our domain has a transport term that we can express thanks to the flow Φ −t (x, v), and the transported datum is f B or f 0 depending on the case the backward trajectory touches a reservoir in the time interval [0, t] on not. We use the function τ to distinguish these two cases. We consider the collision operator 2λη ε K as the source term for the linear problem 9. So we construct the following expression for g ε The notation χ represents the characteristic function.
The meaning of 10 is clear: we separate the contribution given to g ε from trajectories transporting the initial datum f 0 , having no collisions with scatterers and never hitting a reservoirs in the time interval [0, t]; the contribution from trajectories transporting the initial datum f B exiting from a reservoirs at time τ and than moving in Ω until the time t without colliding any scatterers; finally the last term is the contribution due to trajectories having the last collision with a scatterers at time t 1 and never touching the reservoirs in the time interval [t 1 , t].
We iterate the procedure by using 10 again for the g ε in the last integral and from 8 we find: By successive iterations we write the series expansion for the density of particles g ε (x, v, t) as Note that the series are clearly convergent in L +∞ . In 12 the terms with χ(τ = 0) define g in ε that represents the contributions to g ε due to trajectories that stay in Ω for every time in [0, t] while the terms with χ(τ > 0) define g out ε that collects the contributions due to trajectories leaving a mass reservoir at time τ > 0.
Note that g out ε solves the problem 1-3 with initial datum f 0 = 0. We will use the shorthand notation Φ −s (x, v) instead of Φ −s m (x, v, t 1 , . . . , t m , δ 1 , . . . , δ m ) where it is clear by the context to which sequence of collisions we refer. Moreover, the terms with zero collision will be included in the series as the m = 0 terms.
We denote by S ε acting on any h ∈ L ∞ (Ω×S 1 ) the Markov semigroup associated to the g in ε term in 12 for an initial datum h, namely so that in 12 g in ε (t) = S ε (t)f 0 . Proposition 1. There exists ε 0 > 0 such that for any ε < ε 0 and for any h ∈ L ∞ (Ω × S 1 ) it holds Note that in the estimate 14 we are considering t = η ε and the estimate is saying that there is a strictly positive probability for a backward trajectory to exit from Ω in a time of the order of η ε .
Proof of Theorem 2.1. From 12 the stationary solution g S ε of the problem 5 verifies for every t 0 > 0. We can formally write it by iterating the previous one in the form of the Neumann series In order to verify the existence and uniqueness of g S ε we show that the 15 converges. Indeed from Proposition 1 and 15, chosen As a consequence the Neumann series 15 converges in L ∞ and identifies a single element in L ∞ . Choosing an arbitrary t 0 bigger than η ε of the same order of η ε and thanks to the semigroup property of S ε it follows that g S ε does not depend on the time t 0 . So there exists a unique stationary solution g S ε ∈ L ∞ (Ω × S 1 ) satisfying 5.
In order to prove Theorem 2.2 we need some properties of the linear Boltzmann operator L defined in 2. We summarize them in the next lemma.
Lemma 5.1. Let L be the operator defined in 2, then L is a selfadjoint operator on L 2 (S 1 ) and has the form L = 2λ(K − I) where K is a selfadjoint and compact operator (in L 2 (S 1 )). Moreover, K is positive and its spectrum is contained in [0, 1]. The value 0 is the only accumulation point for the spectrum and 1 is a simple eigenvalue. So it holds that {KerL} ⊥ = {h ∈ L 2 (S 1 ) : S 1 dv h(v) = 0} and there exists C > 0 such that for any h ∈ L ∞ (S 1 ) that verifies S 1 dv h(v) = 0 we have Proof of Lemma 5.1. The existence and the estimate of norm of L −1 are discussed in Lemma 4.1 from Section 4.1 of [3]. The compactness of the operator K and the spectral property of L are discussed in [18].
Proof of Theorem 2.2. The proof makes use of the Hilbert expansion (see e.g. [3,6,18]). Assume that g S ε has the following form where g (k) are not depending on η ε . We require g (0) to satisfy the same Dirichlet boundary conditions as the whole solution g S ε on ∂Ω L ∪ ∂Ω R : By imposing that g S ε solves 5 and by comparing terms of the same order we get the following chain of equations: where we used that Lg (0) (x) = 0 since g (0) is independent of v. The first two equations read is an odd function of v, it belongs to (KerL) ⊥ . So we can invert the operator L and set where ζ (1) (x) ∈ Ker L and L −1 (v · ∇ x g 0 ) is an odd function of v since L −1 preserves the parity, namely it maps odd (even) function of v in odd (even) functions (see [18]). We integrate equation (ii) with respect to the uniform measure on S 1 . We can notice that S 1 dv v · ∇ x ζ (1) (x) = 0 (ζ (1) depends only on x, so the function in the integral is odd in the velocity) and S 1 dv Lg (2) = 0 (since operator L preserves mass), so by 18 we obtain 1 2π By expanding the scalar product and using the linearity of L −1 we get We define the 2×2 matrix D i,j = 1 2π S 1 dv v i (−L −1 )v j and we observe that D ij = 0 if i = j as follows by the change v i → −v i while D 11 = D 22 = D > 0 thanks to the isotropy and the spectral property of the operator (see [18]). Hence D is given by the formula 33 and the integrated equation (ii) becomes 1 2π We require g S ε (x, v) to satisfy the reflective boundary condition g S ε (x, v ) = g S ε (x, v) on ∂Ω E . By imposing it on the first term g (1) (x, v) = g (1) (x, v ) for every x ∈ ∂Ω E , v · n < 0, from 18 we obtain By means of the elastic collision rule v = v − 2(v · n)n, the linearity of L −1 allow us to write Left and right members in 22 are the same if and only if (n · ∇ x g (0) )L (−1) (v · n) = 0. Since S 1 dv v·n = 0 we get L (−1) (v·n) = 0, so the only possibility is (n·∇ x g (0) ) = 0. Therefore g (0) (x) has to satisfy the Neumann boundary conditions ∂ n g (0) (x) = 0, for all x ∈ ∂Ω E . From the previous one, 21 and 17 we have shown that the term We can deal with this mixed problem following the method of [28], Chapt. II. Furthermore, regularity results guarantee g 0 ∈ C ∞ (Ω) (see [19], Chapt. 6). Since 19 shows that S 1 dv v · ∇ x g (1) = 0, we can invert L in equation (ii) to obtain (24) where ζ (2) belongs to the kernel of L. Now, integrating the third equation v · ∇ x g (2) (x) = Lg (3) (x, v) with respect to the uniform measure on S 1 , we find thanks to 24 The last integral is null because of the independence of ζ (2) (x) from v. The first integral is null because the function in the integral is an odd function of the velocity thanks to the fact that the operator L −1 preserves the parity. The 25 becomes Since there are no restriction on the choice of the boundary condition, we impose the Dirichlet data ζ (1) (x) = 0 on the boundary ∂Ω L ∪ ∂Ω R . So that by the previous and 26 we find ζ (1) (x) ≡ 0 and hence g (1) Because of the 21 the first term of the right hand side of equation 24 is null too. So 24 reduces to g (2) Moreover from the third equation we get, by inverting L, with ζ (3) (x) belonging to KerL.
By integrating on S 1 the fourth equation v · ∇ x g (3) = Lg (4) and by exploiting that S 1 dv Lg 4 (x, v) = 0 and that We choose zero boundary condition at the reservoirs, namely ζ (2) We can now write the expansion for g S ε as The remainder R ηε satisfies We required g (0) to satisfy the same boundary conditions as the whole solution at contact with the reservoirs, namely on ∂Ω L ∪ ∂Ω R , so the boundary conditions for R ηε read Note that the problem 29-30 has the form of 5. From Theorem 2.1 we know that it admits a unique solution in L ∞ .
From the 28, thanks to the the fact that both g (1) and R ηε are bounded in L ∞ norm, we conclude that g S ε → g (0) .
In order to prove Proposition 1 we follow the strategy of the proof of Proposition 3.1 in [3]. Here we have the additional difficulty of the specular reflective boundaries of horizontal sides of the strip and the presence of the obstacles in Ω. In the proof are exploited the diffusive limit of the linear Boltzmann equation in a L ∞ setting and in a bigger domain containing Ω as stated in Proposition 2 below and the properties of L summarized in Lemma 5.1.
We construct the extended domain Λ as the infinite strip constructed by removing the left and right sides of Ω and keeping the upper and lower elastic boundaries at x 2 = 0 and x 2 = L 2 and the obstacles into Ω (see Figure 15). We call ∂Λ E the union of upper and lower sides of Λ with the obstacles boundaries. We introduce h ε : Λ × S 1 × [0, T ] → R + the solution of the following rescaled linear Boltzmann equation where ρ 0 (x) is a smooth function of the only variable x (local equilibrium).

Proposition 2.
Let h ε be the solution of 31, with an initial datum ρ 0 ∈ C ∞ (Λ) such that there exists M > 0 with ρ 0 (x) = 0 if |x| > M and ∂ n ρ 0 (x) = 0 for x ∈ ∂Λ E . Then, as ε → 0, h ε converges to the solution of the heat equation where the diffusion coefficient D is given by the formula The convergence is in L ∞ ([0, T ]; L ∞ (Λ × S 1 )).
Proof of Proposition 1. The semigroup S ε defined in 13 can be equivalently written as extended to functions belonging to L ∞ (Λ × S 1 ), namely for any f ∈ L ∞ (Λ × S 1 ), where χ Ω is the characteristic function of Ω and Φ −t (x) is the first component (the position) of Φ −t (x, v), the backward flux individuated by x, v, t 1 , . . ., t m , δ 1 , . . ., δ m . The addition of χ Ω (Φ −t (x)) guarantees together with χ(τ = 0) that the dynamics stay internal to Ω. Moreover, the following estimate holds We construct χ δ Ω ∈ C ∞ (Λ), a mollified version of χ Ω , χ δ Ω ≥ χ Ω , χ δ Ω ≤ 1 and Note that the series in 35 defines a function F which solves Defining G ε (x, v, t) as F (x, v, η ε t), G ε solves 31 with initial datum ρ 0 = χ δ Ω . Thanks to Proposition 2 we know that at time t = 1 where ρ δ solves 32 with initial datum χ δ Ω and ω(ε) denotes a positive function vanishing with ε. Moreover, we can notice that the function ρ δ is the solution of a diffusion equation with initial datum 0 ≤ χ δ Ω ≤ 1 with support in a bounded subset of the infinite strip Λ. By the strong maximum principle we know that for the positive time t = 1, it holds that ρ δ (x, 1) < 1. Therefore for ε small enough where we have used 35 for t = η ε .
Proof of Proposition 2. Let h ε : Λ × S 1 × [0, T ] the solution of 31. We use the Hilbert expansion technique to prove that h ε converges to the solution of the heat equation 32. We search h ε of the form with coefficient h (k) not depending on η ε . By imposing that h ε solves 31 and comparing terms of the same order we find the identity Lh (0) (x, t) = 0 and the chain of equations We impose that h (0) satisfy the same initial condition of the whole solution h ε , namely Let us start from the first equation (i) v · ∇ x h (0) = Lh (1) . Thanks to the Fredholm alternative and by proceeding as in the proof of Theorem 2.2, we can solve equation (i) if and only if the left hand side belongs to Since v · ∇ x h (0) (x) is an odd function of v, it belongs to (KerL) ⊥ . So we can invert the operator L finding h (1) (x, v, t) = L −1 (v · ∇ x h (0) (x, t)) + ζ (1) (x, t).
where ζ (1) (x, t) is a function to be determined in the kernel of L. Recall that L −1 preserves the parity.
We integrate the second equation (ii) ∂ t h (0) + v · ∇ x h (1) = Lh (2) with respect to the uniform measure on the sphere S 1 . Thanks to the equation 38 and the observations that S 1 dv Lh (2) = 0 and S 1 dv v · ∇ x ζ (1) (x, t) = 0, it holds As in the proof of Theorem 2.2 defining D i,j = 1 2π S 1 dv v i (−L −1 )v j , we find that the diffusion coefficient D is given by the formula 33 h ε (x, v) has to satisfy the reflective boundary condition h ε (x, v , t) = h ε (x, v, t) on ∂Λ E . By imposing it on the first term h (1) (x, v, t) = h (1) (x, v , t) for every x ∈ ∂Ω E , v · n < 0, we obtain proceeding in the same way of the proof of Theorem 2.2 that h (0) (x, t) has to satisfy the Neumann boundary conditions ∂ n h (0) (x, t) = 0, for all x ∈ ∂Λ E .
We have so shown that the term h (0) (x, t) solves the problem In particular h (0) (t) ∈ L ∞ (Λ × S 1 ) for any t ≥ 0. The equation 40 allow us to verify that when integrating the equation (ii) the left hand side vanishes. It implies that we can invert operator L finding where ζ (2) (x, t) is a function in Ker L.
Observe now that the left hand side of equation (iii) has null integral on S 1 due to 43. By inverting L we obtain the formula for h (3) where ζ (3) ∈ Ker L. We integrate now the equation (iv) ∂ t h (2) + v · ∇ x h (3) = Lh (4) with respect to the uniform measure on S 1 . We find the equation for ζ (2) where the source S(x, t) is given by 1 (x, v, t))).
By means of the series expansion found in 35, the solution can be written in the following way:   T ηε (t) ∞ ≤ C < +∞.
So the remainder is uniformly bounded too. Hence from the estimates and 46 it follows that h ε converges to h (0) in L ∞ for η ε → ∞.