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

Kinetic and Related Models c American Institute of Mathematical Sciences Volume 11, Number 6, December 2018 doi:10.3934/krm.2018058 pp. 1475–1501 LINEAR BOLTZMANN DYNAMICS IN A STRIP WITH LARGE REFLECTIVE OBSTACLES: STATIONARY STATE AND RESIDENCE TIME Alessandro Ciallella∗ and Emilio N. M. Cirillo Dipartimento di Scienze di Base e Applicate per l’Ingegneria Sapienza Università di Roma via A. Scarpa 16, I – 00161, Roma, Italy (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 2010 Mathematics Subject Classification. Primary: 82B21, 76P05, 82B40; Secondary: 60J60. Key words and phrases. Residence time, linear Boltzmann, Kinetic Theory, Monte Carlo methods, obstacles, diffusion. ∗ Corresponding author: Alessandro Ciallella. 1475 1476 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO 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 LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1477 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, L1 ) × (0, L2 ) ⊂ R2 . This strip has two open boundaries, that we think as the left side ∂ΩL = {0} × (0, L2 ) and the right side ∂ΩR = {L1 } × (0, L2 ). 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, L1 ) × {L2 } and lower side (0, L1 ) × {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. Ω ρR ρL ∂ΩL ∂ΩR ∂ΩE Powered by TCPDF (www.tcpdf.org) Figure 1. Domain Ω: strip with large fixed obstacles, where ∂ΩL and ∂ΩR are the vertical open boundaries and ∂ΩE are reflective boundaries. 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 ∈ R2 : |v| = 1}. We will consider the operator L with elastic collision kernel. The equation for (x, v) ∈ Ω×S 1 and positive times t reads Powered by TCPDF (www.tcpdf.org) (∂t + v · ∇x )g(x, v, t) = ηε Lg(x, v, t), x ∈ Ω, v ∈ S 1 , t ≥ 0 (1) 1478 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO where, by the elastic collision rule v 0 = v − 2(n · v)n, the operator L is defined for any f ∈ L1 (S 1 ) as Z 1 [f (v 0 ) − f (v)] dδ. (2) Lf (v) = λ −1 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. v~′ ~ n α α 1 ~v δ Powered by TCPDF (www.tcpdf.org) 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: ( g(x, v 0 , t) = g(x, v, t) x ∈ ∂ΩE , v · n < 0, t ≥ 0 (3) g(x, v, t) = fB (x, v) x ∈ ∂ΩL ∪ ∂ΩR , v · n > 0, t ≥ 0, Powered by TCPDF (www.tcpdf.org) where fB is defined in 4 below and v 0 is given by the elastic collision rule v 0 = 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 f0 (x, v) ∈ L∞ (Ω × S 1 ) and we define fB (not depending on t) as ( ρL /2π x ∈ ∂ΩL , v · n > 0 fB (x, v) := (4) ρR /2π x ∈ ∂ΩR , v · n > 0, LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1479 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:  S S  x ∈ Ω, v ∈ S 1 v · ∇x g (x, v) = ηε Lg S 0 S (5) g (x, v ) = g (x, v) x ∈ ∂ΩE , v · n < 0   S g (x, v) = fB (x, v) x ∈ ∂ΩL ∪ ∂ΩR , v · n > 0. 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  ∆ρ(x) = 0 x∈Ω    ρ(x) = ρ x ∈ ∂ΩL L (6)  ρ(x) = ρR x ∈ ∂ΩR    ∂n ρ(x) = 0 x ∈ ∂ΩE . Theorem 2.1. If ε > 0 is sufficiently small there exists a unique stationary solution gεS ∈ L∞ (Ω × S 1 ) of 5. Theorem 2.2. The stationary solution gεS of 5 verifies gεS → ρ (7) as ε → 0, where ρ(x) is the solution to the problem 6. The convergence is in L∞ (Ω × S 1 ). 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 tm 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 1480 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO 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 th , the time until the hit with a scatterers, following the exponential law of mean tm , with tm a fixed parameter. The particle travels with velocity v for a time th . If during this time it hits the elastic boundary ∂ΩE , its velocity changes performing an elastic collision. At time th we simulate an hit with a scatterers by picking an impact parameter δ uniformly in [−1, 1] and changing the velocity from v to v 0 = 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 f0 , 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 tm 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 tm 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 htm (x). So we fix a cell in contact with the reservoir where we calculated a sojourn time tc and consider the value fB of the stationary solution. We choose as multiplication constant c = fB (x)/tc . So the simulated solution htm (x) is constructed by multiplying the sojourn time in each cell for this constant c. In the sequel we will show that for tm sufficiently small and N big enough, the simulated solution htm (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 · 107 particles. LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1481 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 ρ(x1 , x2 ) = 1 − x1 /8, where we are denoting by (x1 , x2 ) ∈ 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 tm , to understand which values of the mean time tm provide a good approximation of the solution we are looking for. In Figures 3 and 4 we show that the choice of tm of the order of 10−2 is suitable for our purpose. Indeed in Figures 3 we compare the simulated solution htm with the analytic solution for the values tm = 2 · 10−1 , tm = 10−1 , tm = 2 · 10−2 . We see in a 3D plot and in a 2D plot, obtained from the previous one by averaging on the x2 variable, that htm becomes closer to ρ(x1 ) when tm decreases. So in Figures 4 we fix the parameter tm = 10−2 and we verify that htm is close to the function ρ(x1 ) showing the relative error |htm − ρ|/ρ. 1 stationary density stationary density 0.9 1 0.9 0.8 0.7 0.6 0.5 1 0 1 0.5 x1 2 3 0.8 0.7 0.6 0.5 x2 0 4 0 1 2 3 4 x1 Figure 3. Plot of the simulated solutions htm in a 3D plot and in a 2D plot constructed by averaging on the x2 variable: in dark gray tm = 2 · 10−1 , in gray tm = 10−1 , in light gray tm = 2 · 10−2 . In black (grid and dashed line) the analytic solution ρ of the associated Laplace problem. 1 0.016 0.012 x2 0.01 0.5 0.008 0.006 0.004 relative error 0.014 0.002 0 0 0 1 2 3 4 x1 Figure 4. Simulation parameter tm = 10−2 : relative error |htm − ρ|/ρ. 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 tm = 10−2 , since we have shown that in the empty case this choice for the exponential clock allows to 1482 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO construct a numerical solution htm that is close to the analytical solution ρ(x1 , x2 ) 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 ρ(x1 , x2 ) of the associated Laplace problem 6. We show the plots of the htm and ρ and the map of the relative error |htm − ρ|/ρ 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. 0.014 0.01 0.008 0.5 0.006 0.004 1 0 1 0.5 x1 2 3 0.002 0 x2 0 0 1 4 0 2 3 4 x1 1 0.012 0.008 0.5 0.006 0.004 relative error 0.01 1 0.9 0.8 0.7 0.6 0.5 x2 stationary density relative error 0.012 x2 stationary density 1 1 0.9 0.8 0.7 0.6 0.5 0.002 1 0 1 0.5 x1 2 3 4 0 x2 0 0 0 1 2 3 4 x1 Figure 5. Simulation parameter tm = 10−2 : on the left in gray the numerical solution htm and in black the solution ρ of the associated Laplace problem; on the right the relative error |htm − ρ|/ρ. Into the strip there is a square obstacle with side 8 · 10−1 . 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 x2 anymore as it was in the empty strip case. However, we can notice that before and beyond the obstacles in the x1 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 {x1 } × (0, 1) ∩ Ω does not depend on x1 . Indeed, this should follow from the Fick’s law, that we expect to be valid also in presence of obstacles (in absence LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1483 0.014 0.01 0.008 0.5 0.006 0.004 0.002 1 0 1 0.5 x1 2 3 x2 relative error 0.012 x2 stationary density 1 1 0.9 0.8 0.7 0.6 0.5 0 0 0 4 0 1 2 3 4 x1 Figure 6. Simulation parameter tm = 10−2 : on the left in gray the numerical solution htm and in black the solution ρ of the associated Laplace problem; on the right the relative error |htm − ρ|/ρ. In the strip is placed a very thin obstacle with height of 0.8. 0.016 0.012 0.01 0.5 0.008 0.006 0.004 1 1 0.5 x1 2 3 0.002 0 x2 0 0 1 4 0 2 3 4 x1 x2 stationary density 1 1 0.9 0.8 0.7 0.6 0.5 1 0 1 0.5 x1 2 3 4 0 x2 0.018 0.016 0.014 0.012 0.01 0.008 0.006 0.004 0.002 0 0.5 0 0 1 2 3 relative error 0 relative error 0.014 x2 stationary density 1 1 0.9 0.8 0.7 0.6 0.5 4 x1 Figure 7. Simulation parameter tm = 10−2 : on the left in gray the numerical solution htm and in black the solution ρ of the associated Laplace problem; on the right the relative error |htm − ρ|/ρ. In the first line we show the case of two squared obstacles with side 6 · 10−1 , in the second one a couple of rectangular obstacles, taller and thinner than the squares. 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 {x1 = 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. 1484 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO 4. Residence time. We consider the domain (0, L1 )×(0, L2 ) 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 tm 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). LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 420 400 395 410 390 400 residence time residence time 1485 390 380 370 385 380 375 370 365 360 360 355 350 350 0 0.2 0.4 0.6 0.8 1 0 0.2 obstacle height 0.6 0.8 1 Ω Ω Powered by TCPDF (www.tcpdf.org) 0.4 obstacle height Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) 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: L1 = 4, L2 = 1, tm = 2 · 10−2 , total number of inserted particles 108 , the total number of particles exiting through the right boundary varies from 5.3 · 105 to 3.6 · 105 (on the left) and from 5.3 · 105 to 2.1 · 105 (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. ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO 365 370 360 360 residence time residence time 1486 355 350 345 340 350 340 330 320 310 335 300 0 0.2 0.4 0.6 0.8 1 0 0.2 obstacle height Powered by TCPDF (www.tcpdf.org) 0.6 0.8 1 Ω Ω Powered by TCPDF (www.tcpdf.org) 0.4 obstacle height Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) 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: L1 = 4, L2 = 1, tm = 2 · 10−2 , total number of inserted particles 108 , the total number of particles exiting through the right boundary varies from 5.2 · 105 to 1.4 · 105 (on the left) and from 5.2 · 105 to 1.1 · 105 (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, L1 ) × (0, L2 ) 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, L1 ) × (0, L2 ) \ (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, 400 365 380 360 360 355 residence time residence time LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 340 320 300 280 1487 350 345 340 335 260 330 0 0.8 1.6 2.4 3.2 4 0 0.2 obstacle width 0.4 0.6 0.8 1 obstacle side Ω Ω Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Figure 10. Residence time vs. width of a centered rectangular obstacle with fixed height 0.8 (on the left) and vs. the side length of a centered squared obstacle (on the right). Simulation parameters: L1 = 4, L2 = 1, tm = 2·10−2 , total number of inserted particles 108 , the total number of particles exiting through the right boundary varies from 4.2 · 105 to 1.1 · 105 (on the left) and from 5.3 · 105 to 1.3 · 105 (on the right) depending on the obstacle width. The solid lines represent the value of the residence time measured for the empty strip (no obstacle). Powered by TCPDF (www.tcpdf.org) 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 1488 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO 500 480 480 460 residence time residence time 460 440 420 400 380 360 440 420 400 380 340 360 320 0.7 1.4 2.1 2.8 3.5 0.7 obstacle position Ω C 1.4 2.1 2.8 3.5 obstacle position Ω C Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Figure 11. Residence time vs. position of the center of the obstacle. The obstacle is a square of side length 0.8 on the left and a rectangle of side lengths 0.04 and 0.8 on the right. Simulation parameters: L1 = 4, L2 = 1, tm = 2 · 10−2 , total number of inserted particles 108 , the total number of particles exiting through the right boundary is stable at the order of 2.6 · 105 (on the left) and of 4 · 105 (on the right) not depending on the obstacle position. The solid lines represent the value of the residence time measured for the empty strip (no obstacle). 5 4 3 2 1 0 200 residence time local residence time Powered by TCPDF (www.tcpdf.org) 120 80 40 1 x2 160 0.5 0 0 0 2 1 x1 3 4 0 0.2 0.4 0.6 0.8 obstacle height Figure 12. As in the right panel in Figure 8. In the left panel the height of the obstacle is equal to 0.8. Left panel: the mean time spent by particles crossing the strip in each point of the strip (0.02× 0.02 cells have been considered) for the empty strip case (black) and in presence of the obstacle (gray). Right panel: residence time in regions L (circles), C (squares), and R (triangles) in presence of the obstacle (gray) and for the empty strip case (black). 1 300 240 180 120 60 1 x2 1489 360 10 8 6 4 2 0 residence time local residence time LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 0 0.5 0 0 1 2 3 0 4 1 2 3 4 obstacle width x1 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. 400 4 residence time local residence time 500 6 2 0 200 100 1 x2 300 0.5 0 0 0 1 2 3 4 0.8 x1 1.6 2.4 3.2 obstacle position Figure 14. As in 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) = f0 (x, v) ∈ L∞ (Ω × 1 S ). We can express the operator defined in 2 as Lf (v) = 2λ(K − I)f (v) , where Z 1 1 dδ f (v 0 ) (8) (Kf )(v) = 2 −1 and I is the identity. Therefore the equation 1 can be written as ∂t f + (v · ∇x + 2ηε λI)f = 2ηε λKf. (9) 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 R2 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 e−t v·∇ f (x, v) = f (x − vt, v). We want to consider the semigroup generated by A on our domain Ω initial datum f0 and boundary conditions 3. Recall that ∂ΩE is a specular reflective 1490 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO boundary while on ∂ΩL ∪ ∂ΩR the system is in contact with reservoirs with particle densities fB (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 ti , δ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 − t1 it does not hit any scatterers, but if the particle reaches the elastic boundary ∂ΩE during this time, the velocity −v becomes −v 0 following the elastic collision rule −v 0 = −(v − 2(n · v)n), where n is the inward pointing normal to Ω. After a time t−t1 the particle performs a collision with impact parameter δ1 that produces the velocity −v1 . Then again the particles travels for a time t1 − t2 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, t1 , . . ., tm , δ1 , . . ., δm , we define the flow Φ−t m (x, v, t1 , . . . , tm , δ1 , . . . , δm ) as the backward trajectory starting from x with velocity v and having m transition in velocity obtained after a time t − t1 , . . ., ti − ti+1 , . . ., tm (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 0 = v − 2(n · v)n any time the elastic boundary ∂ΩE is reached. We define the function τ = τ (x, v, t, t1 , . . . , tm , δ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 t1 , . . . , tm , δ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 fB or f0 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ε gε (x, v, t) = e−2ληε t f0 (Φ−t 0 (x, v))χ(τ = 0) −(t−τ ) (x, v))χ(τ > 0) + e−2ληε (t−τ ) fB (Φ0 Z t −(t−t1 ) (x, v), t1 )χ(τ < t1 ) dt1 . + ληε e−2ληε (t−t1 ) 2Kgε (Φ0 (10) 0 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 f0 , having no collisions with scatterers and never hitting a reservoirs in the time interval [0, t]; the contribution from trajectories transporting the initial datum fB 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 t1 and never touching the reservoirs in the time interval [t1 , t]. We iterate the procedure by using 10 again for the gε in the last integral and from 8 we find: LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1491 −(t−τ ) −2ληε (t−τ ) gε (x, v, t) = e−2ληε t f0 (Φ−t fB (Φ0 (x, v))χ(τ > 0) 0 (x, v))χ(τ = 0) + e  Z 1 Z t dδ1 e−2ληε t1 f0 (Φ−t + ληε dt1 e−2ληε (t−t1 ) 1 (x, v, t1 , δ1 ))χ(τ = 0) 0 +e −2ληε (t1 −τ ) t1 Z −1 −(t−τ ) fB (Φ1 (x, v, t1 , δ1 ))χ(τ > 0)χ(τ < t1 ) dt2 e−2ληε (t1 −t2 ) + ληε 0 Z 1 × −1 −(t−t2 ) dδ2 gε (Φ1  (x, v, t1 , δ1 ), t2 )χ(τ < t2 ) . (11) By successive iterations we write the series expansion for the density of particles gε (x, v, t) as gε (x, v, t) = e−2ληε t f0 (Φ−t 0 (x, v))χ(τ = 0) Z t Z X −2ληε t m + e (ληε ) dt1 . . . 0 m≥1 Z 1 Z 1 Z 1 tm−1 dtm 0 dδm χ(τ = 0)f0 (Φ−t m (x, v, t1 , . . . , tm , δ1 , . . . , δm )) −1 −(t−τ ) e−2ληε (t−τ ) fB (Φ0 (x, v))χ(τ > 0) Z Z t tm−1 X (ληε )m dt1 . . . dtm 0 0 m≥1 × dδ1 . . . −1 + + Z 1 × dδm χ(τ < tm )χ(τ > 0)e−2ληε (t−τ ) dδ1 . . . −1 −1 ) × fB (Φ−(t−τ (x, v, t1 , . . . , tm , δ1 , . . . , δm )) m = gεin (x, v, t) + gεout (x, v, t). (12) +∞ 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 f0 = 0. We will use the shorthand notation Φ−s (x, v) instead of Φ−s m (x, v, t1 , . . . , tm , δ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 Z t Z tm−1 X (Sε (t)h)(x, v) = e−2ληε t (ληε )m dt1 . . . dtm 0 m≥0 Z 1 × (13) 1 dδ1 . . . −1 so that in 12 gεin (t) = Sε (t)f0 . 0 Z −t dδm χ(τ = 0)h(Φ (x, v)), −1 1492 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO Proposition 1. There exists ε0 > 0 such that for any ε < ε0 and for any h ∈ L∞ (Ω × S 1 ) it holds kSε (ηε )hk∞ ≤ αkhk∞ , α < 1. (14) 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 gεS = gεout (t0 ) + Sε (t0 )gεS , for every t0 > 0. We can formally write it by iterating the previous one in the form of the Neumann series X gεS = (Sε (t0 ))N gεout (t0 ). (15) N ≥0 In order to verify the existence and uniqueness of gεS we show that the 15 converges. Indeed from Proposition 1 and 15, chosen t0 = ηε X X kgεS k∞ ≤ k(Sε (ηε ))N gεout k∞ ≤ αN kgεout (ηε )k∞ N ≥0 N >0 1 1 kg out k∞ ≤ max{ρL , ρR }. ≤ 1−α ε 1−α As a consequence the Neumann series 15 converges in L∞ and identifies a single element in L∞ . Choosing an arbitrary t0 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 t0 . 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 L2 (S 1 ) and has the form L = 2λ(K − I) where K is a selfadjoint and compact operator (in L2 (S 1 )). Moreover, K is positive and its spectrum is contained in [0, 1]. The value 0 is the only accumulation point for theRspectrum and 1 is a simple eigenvalue. So it holds that {KerL}⊥ = {h ∈ L2 (S 1 ) : RS 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 kL−1 hk∞ ≤ Ckhk∞ . (16) 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 k +∞  X 1 gεS (x, v) = g (0) (x) + g (k) (x, v), ηε k=1 LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1493 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 : ( g (0) (x) = ρL x ∈ ∂ΩL (17) (0) g (x) = ρR x ∈ ∂ΩR . By imposing that gεS solves 5 and by comparing terms of the same order we get the following chain of equations: v · ∇x g (k) = Lg (k+1) , k ≥ 0, where we used that Lg (0) (x) = 0 since g (0) is independent of v. The first two equations read (i) v · ∇x g (0) (x) = Lg (1) (x, v), (ii) v · ∇x g (1) (x, v) = Lg (2) (x, v). Let us consider the first one. By the Fredholm alternative, this equation has a solution if and only if the left hand side belongs to (KerL)⊥ . We recall that the null space of L is constituted by the constant functions (with respect to v), so we can solve equationR(i) if and only if the left hand side belongs to (KerL)⊥ = {h ∈ L2 (S 1 ) such that S 1 dv h(v) = 0} (see Lemma 5.1). Since v · ∇x g (0) (x) is an odd function of v, it belongs to (KerL)⊥ . So we can invert the operator L and set g (1) (x, v) = L−1 (v · ∇x g (0) (x)) + ζ (1) (x), (18) 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 R notice that S 1 dv v · ∇x ζ (1) (x) = 0 (ζ (1) depends only on x, so the function in the R integral is odd in the velocity) and S 1 dv Lg (2) = 0 (since operator L preserves mass), so by 18 we obtain Z  1 dv v · ∇x (L−1 (v · ∇x g (0) (x))) = 0. (19) 2π S1 By expanding the scalar product and using the linearity of L−1 we get − 2 X Di,j ∂xi ∂xj g (0) (x) = 0. (20) i,j=1 R 1 dv vi (−L−1 )vj and we observe that Dij = 0 We define the 2×2 matrix Di,j = 2π S1 if i 6= j as follows by the change vi → −vi while D11 = D22 = D > 0 thanks to the isotropy and the spectral property of the operator (see [18]). Hence D is given by the formula 33 Z 1 D= dv v · (−L)−1 v, 4π S 1 and the integrated equation (ii) becomes Z  1 dv v · ∇x (L−1 (v · ∇x g (0) (x))) = 0 ⇔ −D∆x g (0) (x) = 0. (21) 2π S1 We require gεS (x, v) to satisfy the reflective boundary condition gεS (x, v 0 ) = on ∂ΩE . By imposing it on the first term g (1) (x, v) = g (1) (x, v 0 ) for every gεS (x, v) 1494 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO x ∈ ∂ΩE , v · n < 0, from 18 we obtain L(−1) (v · ∇x g (0) ) + ζ (1) (x) = L(−1) (v 0 · ∇x g (0) ) + ζ (1) (x). (22) By means of the elastic collision rule v 0 = v − 2(v · n)n, the linearity of L−1 allow us to write L(−1) ((v − 2(v · n)n) · ∇x g (0) ) = L(−1) (v · ∇x g (0) ) − 2(n · ∇x g (0) )L(−1) (v · n). Left and are the same if and only if (n · ∇x g (0) )L(−1) (v · n) = 0. R right members in 22(−1) Since S 1 dv v·n = 0 we get L (v·n) 6= 0, so the only possibility is (n·∇x g (0) ) = 0. (0) Therefore g (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 g (0) (x) solves the problem  ∆x g (0) (x) = 0 x∈Ω    g (0) (x) = ρ x ∈ ∂ΩL L (23) (0)  g (x) = ρ x ∈ ∂ΩR R    ∂n g (0) (x) = 0 x ∈ ∂ΩE . We can deal with this mixed problem following the method of [28], Chapt. II. Furthermore, regularity Rresults 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 g (2) (x, v) = L−1 (v · ∇x L−1 (v · ∇x g (0) (x))) + L−1 (v · ∇x ζ (1) (x)) + ζ (2) (x), (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 Z dv v · ∇x (L−1 (v · ∇x L−1 (v · ∇x g (0) (x))))+ S1 Z Z (25) −1 (1) (2) + dv v · ∇x (L (v · ∇x ζ (x))) + dv v · ∇x (ζ (x)) = 0. S1 S1 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 Z dv v · ∇x (L−1 (v · ∇x ζ (1) (x))) = −D∆x ζ (1) (x) = 0 (26) S1 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) (x, v) = L−1 (v · ∇x g (0) (x)). Because of the 21 the first term of the right hand side of equation 24 is null too. So 24 reduces to g (2) (x, v) = ζ (2) (x). Moreover from the third equation we get, by inverting L, g (3) (x, v) = L−1 (v · ∇x g (2) (x, v)) + ζ (3) (x) = L−1 (v · ∇x ζ (2) (x)) + ζ (3) (x), with ζ (3) (x) belonging to KerL. LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1495 ByR integrating on S 1 the fourthRequation v · ∇x g (3) = Lg (4) and by exploiting that S 1 dv Lg 4 (x, v) = 0 and that S 1 dv v · ∇x ζ (3) (x) = 0 we find Z dv v · ∇x (L−1 (v · ∇x ζ (2) (x))) = −D∆x ζ (2) (x) = 0. (27) S1 We choose zero boundary condition at the reservoirs, namely ζ (2) (x) = 0 on ∂ΩL ∪ ∂ΩR , so we find ζ (2) (x) ≡ 0. Then g (2) (x, v) ≡ 0. We can now write the expansion for gεS as gεS = g (0) + 1 1 (1) g + R ηε . ηε ηε (28) The remainder Rηε satisfies v · ∇x Rηε = ηε LRηε . (29) 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 ( Rηε (x, v) = −L−1 (v · ∇x g (0) (x)) x ∈ ∂ΩL ∪ ∂ΩR , v · n(x) > 0, (30) Rηε (x, v 0 ) = Rηε (x, v) x ∈ ∂ Ω̃E , v · n(x) < 0. 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 x2 = 0 and x2 = L2 and the obstacles into Ω (see Figure 15). We call ∂ΛE the union of upper and lower sides of Λ with the obstacles boundaries. Λ ∂ΛE Powered by TCPDF (www.tcpdf.org) Figure 15. Domain Λ: infinite strip with big fixed obstacles: the whole boundaries of Λ is a specular reflective boundary. Powered by TCPDF (www.tcpdf.org) 1496 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO We introduce hε : Λ × S 1 × [0, T ] → R+ the solution of the following rescaled linear Boltzmann equation  2  x∈Λ (∂t + ηε v · ∇x )hε = ηε Lhε 0 (31) hε (x, v , t) = hε (x, v, t) x ∈ ∂ΛE , v · n < 0, t ≥ 0   hε (x, v, 0) = ρ0 (x) x ∈ Λ, 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   x∈Λ ∂t ρ − D∆ρ = 0 (32) ρ(x, 0) = ρ0 (x) x∈Λ   ∂n ρ(x, t) = 0 x ∈ ∂ΛE , t ≥ 0, where the diffusion coefficient D is given by the formula Z 1 D= dv v · (−L)−1 v. 4π S 1 (33) 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 Z t Z m−1 X −2ληε t m (Sε (t)f )(x, v) =χΩ (x) e (ληε ) dt1 . . . dtm 0 m≥0 Z 1 Z dδ1 . . . −1 0 (34) 1 −t −t dδm χ(τ = 0)f (Φ (x, v))χΩ (Φ (x)), −1 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, t1 , . . ., tm , δ1 , . . ., δm . The addition of χΩ (Φ−t (x)) guarantees together with χ(τ = 0) that the dynamics stay internal to Ω. Moreover, the following estimate holds Z t Z m−1 X Sε (t)f ≤kf k∞ e−2ληε t (ληε )m dt1 . . . dtm 0 m≥0 1 Z × 0 1 Z dδm χΩ (Φ−t (x)). dδ1 . . . −1 −1 We construct χδΩ ∈ C ∞ (Λ), a mollified version of χΩ , χδΩ ≥ χΩ , χδΩ ≤ 1 and Ω ⊂ supp(χδΩ ) ⊂ (−δ, L1 + δ) × [0, L2 ]. So we can write Z t Z m−1 X Sε (t)f ≤kf k∞ e−2ληε t (ληε )m dt1 . . . dtm 0 m≥0 Z 1 × dδ1 . . . −1 0 Z 1 dδm χδΩ (Φ−t (x)). −1 (35) LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES Note that the series in 35 defines a function F which solves   x∈Λ (∂t + v · ∇x )F (x, v, t) = ηε LF (x, v, t) F (x, v 0 , t) = F (x, v, t) x ∈ ∂ΛE , v · n < 0, t ≥ 0   F (x, v, 0) = χδΩ (x) x ∈ Λ. 1497 (36) 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 kGε (1) − ρδ (1)k∞ ≤ ω(ε) 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 kSε (ηε )f k∞ ≤kf k∞ kSε (ηε )χδΩ k∞ ≤ kf k∞ (kGε (1) − ρδ (1)k∞ + kρδ (1)k∞ ) ≤kf k∞ (ω(ε) + kρδ (1)k∞ ) < αkf k∞ , α < 1, (37) 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 k +∞  X 1 hε (x, v, t) = h(0) (x, t) + h(k) (x, v, t), ηε k=1 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 v · ∇x h(0) =Lh(1) ∂t h(k) + v · h(k+1) =Lh(k+2) for k ≥ 0. We impose that h(0) satisfy the same initial condition of the whole solution hε , namely h(0) (x, 0) = ρ0 (x). 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 Z (KerL)⊥ = {h ∈ L2 (S 1 ) such that dv h(v) = 0}. S1 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). (38) where ζ (x, t) is a function to be determined in the kernel of L. Recall that L−1 preserves the parity. (1) 1498 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO We integrate the second equation (ii) ∂t h(0) + v · ∇x h(1) = Lh(2) with respect to the uniform measure on the sphereR S 1 . Thanks to the equation 38 and the R observations that S 1 dv Lh(2) = 0 and S 1 dv v · ∇x ζ (1) (x, t) = 0, it holds Z 1 dv ∂t h(0) (x, t) + v · ∇x (L−1 v · ∇x h(0) (x, t)) = 0. (39) 2π S 1 R 1 As in the proof of Theorem 2.2 defining Di,j = 2π dv vi (−L−1 )vj , we find that S1 the diffusion coefficient D is given by the formula 33 Z 1 dv v · (−L)−1 v D= 4π S 1 so that the heat equation for h(0) is ∂t h(0) − D∆x h(0) = 0. (40) 0 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 0 , 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  (0) (0)  x∈Λ ∂t h − ∆x h = 0 (0) (41) h (x, 0) = ρ0 (x) x∈Λ   ∂n h(0) (x, t) = 0 x ∈ ∂ΛE . 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 h(2) (x, v, t) = L−1 (∂t h(0) (x, t) + v · ∇x (L−1 (v · ∇x h(0) (x, t))) (42) + v · ∇x ζ (1) (x, t)) + ζ (2) (x, t), where ζ (2) (x, t) is a function in Ker L. Next equation is (iii) ∂t h(1) + v · ∇x h(2) = Lh(3) . When integrating it with respect to the uniform measure on S 1 , we exploit the fact that the operator L−1 preserves the parity. So, substituting h(1) and h(2) with their expressions given by 38 and 42, the only terms surviving give the equation for ζ (1) ∂t ζ (1) (x, t) − D∆x ζ (1) (x, t) = 0. (43) (1) Since there are no restrictions on the choice of the initial condition for ζ , we fix ζ (1) (x, 0) = 0. So ζ (1) (x, t) ≡ 0 for any (x, t) and the expression for h(1) reduces to h(1) (x, v, t) = L−1 (v · ∇x h(0) (x, t)). By the Lemma 5.1 and the smoothness of h(0) we have sup kh(1) (t)k∞ ≤ C sup k∇x h(0) (t)k∞ < +∞. t∈[0,T ] t∈[0,T ] In the same way, by Lemma 5.1 and smoothness of h(0) it follows that the (2) first term in the expression of h(2) , i.e. h1 = L−1 (∂t h(0) (x, t) + v · ∇x (L−1 (v · (0) ∞ ∞ 1 ∇x h (x, t)))), is in L ([0, T ]; L (Λ × S )), as well as its spatial derivatives. LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1499 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) h(3) (x, v, t) = L−1 (∂t h(1) + v · ∇x h(2) (x, v, t)) + ζ (3) (x, t) (2) = L−1 (∂t L−1 (v · ∇x h(0) (x, t)) + v · ∇x (h1 (x, v, t) +ζ (2) (x, t))) + ζ (3) (44) (x, t), 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) (x, t) ∂t ζ (2) − D∆x ζ (2) = S(x, t), (45) where the source S(x, t) is given by Z 1 dv v · ∇x L−1 (∂t L−1 (v · ∇x h(0) (x, t))) S(x, t) = − 2π S 1 Z 1 (2) − dv v · ∇x L−1 (v · ∇x h1 (x, v, t))). 2π S 1 We consider as initial datum ζ (2) (x, 0) = 0, so we have ζ (2) ∈ L∞ ([0, T ]; L∞ (Λ)) and its spatial derivative as well, since S ∈ L∞ ([0, T ]; L∞ (Λ)). We write the the expansion truncated at order ηε−2 for the solution: hε (x, v, t) = h(0) (x, t) + 1 (1) 1 1 h (x, v, t) + 2 h(2) (x, v, t) + Rηε (x, v, t). ηε ηε ηε (46) We have shown that h(i) (t) ∈ L∞ (Λ × S 1 ) for i = 0, 1, 2. Now we have to prove that even,the remainder Rηε is in L∞ . The remainder Rηε satisfies the equation (∂t + ηε v · ∇x )Rηε = ηε2 LRηε − Tηε (47) with initial condition Rηε (x, v, 0) = −h(1) (x, v, 0) − 1 (2) h (x, v, 0) ηε and boundary conditions Rηε (x, v 0 , t) = Rηε (x, v, t) x ∈ ∂ΛE , v · n < 0. The term Tηε on the left hand side of 47 is Tηε = ∂t h(1) + η1ε ∂t h(2) + v · ∇x h(2) . So Tηε ∈ L∞ ([0, T ]; L∞ (Λ × S 1 )) and thanks to the smoothness hypothesis on ρ0 also the initial datum Rηε (x, v, 0) belongs to L∞ . By denoting by Sηε (t) the semigroup associated to the generator ηε (v · ∇x − ηε L) with reflective boundary conditions on ∂ΛE , the equation 47 becomes Z t ds Sηε (t − s)Tηε (s). Rηε (t) = Sηε (t)Rηε (0) + 0 1500 ALESSANDRO CIALLELLA AND EMILIO N. M. CIRILLO By means of the series expansion found in 35, the solution can be written in the following way: Z ηε t Z m−1 X 2 dt1 . . . dtm Rηε (x, v, t) = e−2ληε t (ληε )m 0 m≥0 Z 1 × Z t ds 0 X 1 dδm Rηε (0)(Φ−ηε t (x, v)) dδ1 . . . −1 −1 + 0 Z 2 e−2ληε (t−s) (ληε )m Z 1 Z dδ1 . . . −1 ηε (t−s) Z dt1 . . . 0 m≥0 × Z m−1 dtm 0 1 dδm Tηε (0)(Φ−ηε (t−s) (x, v), s). −1 Therefore we can estimate sup kRηε (t)k∞ ≤ kRηε (0)k∞ + T sup kTηε (t)k∞ ≤ C < +∞. t∈[0,T ] t∈[0,T ] So the remainder is uniformly bounded too. Hence from the estimates and 46 it follows that hε converges to h(0) in L∞ for ηε → ∞. Acknowledgments. We thank Daniele Andreucci for many useful discussions on mixed boundary problems and for having suggested the reference [28]. We are grateful to Mario Pulvirenti for many illuminating discussions on the topic of the paper and for having suggested the idea of studying the residence time problem in the framework of Kinetic Theory. REFERENCES [1] G. Albi, M. Bongini, E. Cristiani and D. Kalise, Invisible control of self-organizing agents leaving unknown environments, SIAM J. Appl. Math., 76 (2016), 1683–1710. [2] F. Alonso-Marroquin, S.I. Azeezullah, S.A. Galindo-Torres and L.M. Olsen-Kettle, Bottlenecks in granular flow: When does an obstacle increase the flow rate in an hourglass?, Phys. Rev. E , 85 (2012), 020301. [3] G. Basile, A. Nota, F. Pezzotti and M. Pulvirenti, Derivation of the Fick’s law for the Lorentz model in a low density regime, Comm. Math. Phys., 336 (2015), 1607–1636. [4] N. Bellomo and C. Dogbe, On the modeling of traffic and crowds: A survey of models, speculations, and perspectives, SIAM Review , 53 (2011), 409–463. [5] D. Braess, A. Nagurney and T. Wakolbinger, On a paradox of traffic planning, Transportation Science, 39 (2005), 446–450. [6] C. Cercignani, R. Illner and M. Pulvirenti, The Mathematical Theory of Dilute Gases, Springer-Verlag, Berlin, 1994. [7] A. Ciallella, On the linear Boltzmann transport equation: A Monte Carlo algorithm for stationary solutions and residence times in presence of obstacles, in AIMETA 2017 - Proceedings of the 23rd Conference of the Italian Association of Theoretical and Applied Mechanics, 5 (2017), 952–960. [8] A. Ciallella, E. N. M. Cirillo and J. Sohier, Residence time of symmetric random walkers in a strip with large reflective obstacles, Phys. Rev. E , 97 (2018), 052116. [9] E. N. M. Cirillo and M. Colangeli, Stationary uphill currents in locally perturbed zero-range processes, Phys. Rev. E , 96 (2017), 052137. [10] E. N. M. Cirillo, M. Colangeli and A. Muntean, Does communication enhance pedestrians transport in the dark?, Comptes Rendus Mecanique, 344 (2016), 19–23. [11] E. N. M. Cirillo, O. Krehel, A. Muntean and R. van Santen, Lattice model of reduced jamming by a barrier, Phys. Rev. E , 94 (2016), 042115. [12] E. N. M. Cirillo, O. Krehel, A. Muntean, R. van Santen and A. Sengar, Residence time estimates for asymmetric simple exclusion dynamics on strips, Phys. A, 442 (2016), 436–457. LINEAR BOLTZMANN DYNAMICS IN PRESENCE OF OBSTACLES 1501 [13] E. N. M. Cirillo and A. Muntean, Can cooperation slow down emergency evacuations? Comptes Rendus Mécanique, 340 (2012), 625–628. [14] E. N. M. Cirillo and A. Muntean, Dynamics of pedestrians in regions with no visibility–a lattice model without exclusion, Phys. A, 392 (2013), 3578–3588. [15] E. Cristiani and D. Peri, Handling obstacles in pedestrian simulations: Models and optimization, Appl. Math. Model., 45 (2017), 285–302. [16] A. J. Ellery, M. J. Simpson, S. W. McCue and R. E. Baker, Characterizing transport through a crowded environment with different obstacle sizes, The Journal of Chemical Physics, 140 (2014), 054108. [17] R. Escobar and A. De La Rosa, Architectural Design for the Survival Optimization of Panicking Fleeing Victims, in Advances in Artificial Life. ECAL 2003 (eds. W. Banzhaf, J. Ziegler, T. Christaller, P. Dittrich, and J.T. Kim), 2801 Springer (2003), 97–106. [18] R. Esposito and M. Pulvirenti, From Particles to Fluids, Hand-Book of Mathematical Fluid Dynamics Vol.III, North-Holland, Amsterdam, (2004), 1–82. [19] L. C. Evans, Partial Differential Equations, Second edition. Graduate Studies in Mathematics, 19. American Mathematical Society, Providence, RI, 2010. [20] B. W. Fitzgerald, J. T. Padding and R. van Santen, Simple diffusion hopping model with convection, Phys. Rev. E , 95 (2017), 013307. [21] D. Helbing, Traffic and related self-driven many-particle systems, Rev. Mod. Phys., 73 (2001), 1067–1141. [22] D. Helbing, L. Buzna, A. Johansson and T. Werner, Self-organized pedestrian crowd dynamics: Experiments, simulations, and design solutions, Transportation Science, 39 (2005), 1–24. [23] D. Helbing, I. Farkas, P. Molnàr and T. Vicsek, Simulation of pedestrian crowds in normal and evacuation situations, in Pedestrian and Evacuation Dynamics (eds. M. Schreckenberg and S. D. Sharma), Springer, (2002), 21–58. [24] D. Helbing, I. J. Farkas and T. Vicsek, Simulating dynamical features of escape panic, Nature, 407 (2000), 487–490. [25] D. Helbing, P. Molnár, I. J. Farkas and K. Bolay, Self-organizing pedestrian movement, Environment and Planning B: Planning and Design, 28 (2001), 361–383. [26] F. Höfling and T. Franosch, Anomalous transport in the crowded world of biological cells, Rep. Progr. Phys., 76 (2013), 046602, 50 pp. [27] R. L. Hughes, The flow of human crowds, Annual Review of Fluid Mechanics, 35 (2003), 169–182. [28] O. A. Ladyzhenskaya, The Boundary Value Problems of Mathematical Physics, SpringerVerlag, Berlin, 1985. [29] M. Matsumoto and T. Nishimura, Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator, ACM Trans. on Modeling and Computer Simulation, 8 (1998), 3–30. [30] M. Matsumoto and T. Nishimura. A Nonempirical Test on the Weight of Pseudorandom Number Generators, in: Monte Carlo and Quasi-Monte Carlo methods 2000 (eds. K.T. Fang, F.J.Hickernel, and H. Niederreiter), Springer-Verlag, (2002), 381–395. [31] M. A. Mourão, J. B. Hakim and S. Schnell, Connecting the dots: The effects of macromolecular crowding on cell physiology, Biophysical Journal, 107 (2017), 2761–2766. [32] M. J. Saxton, Anomalous diffusion due to obstacles: A Monte Carlo study, Biophysical Journal, 66 (1994), 394–401. [33] K. To, P. Y. Lai and H. K. Pak, Jamming of granular flow in a two-dimensional hopper, Phys. Rev. Lett., 86 (2001), 71–74. [34] I. Zuriguel, A. Garcimartı́n, D. Maza, L. A. Pugnaloni and J. M. Pastor, Jamming during the discharge of granular matter from a silo, Phys. Rev. E , 71 (2005), 051303. [35] I. Zuriguel, A. Janda, A. Garcimartı́n, C. Lozano, R. Arévalo and D. Maza, Silo clogging reduction by the presence of an obstacle, Phys. Rev. Lett., 107 (2011), 278001. Received for publication August 2017. E-mail address: alessandro.ciallella@uniroma1.it E-mail address: emilio.cirillo@uniroma1.it