The BGK equation as the limit of an $N$ particle system

The spatially homogeneous BGK equation is obtained as the limit if a model of a many particle system, similar to Mark Kac's charicature of the spatially homogeneous Boltzmann equation.


INTRODUCTION
The BGK-equation is named after P.L. Bhatnagar, E.P. Gross, and M. Krook, who first presented it in an influential paper published in 1954 [3]. In its original form it is Here f = f (x, v, t) gives the number density of particles in phase-space (x, v) ∈ 3 × 3 . The constant σ > 0 controls collision rate of particles, and Φ = Φ q,T is the Maxwellian distribution where the n(x, t), q(x, t), and T (x, t) represent the local number density, mean velocity and temperature respectivley: The same kind of equation was formulated indpendently by Welander [23]. In [3], one considers charged particles, and E is the electric field computed from the particle density. It is a model of the kinetic Boltzmann equation with the purpose of providing a numerically tractable model, while retaining the most important aspects of the original Boltzmann equation: conservation of mass, momentum and energy, convergence to a unique equilibrium state, monotonicity of entropy, etc. And while easier from a computational point of view, it is considerably more difficult to analyse mathematically, and most theoretical results concerning existence and uniqueness of solutions to the BGK eqution actually hold for a modified version where the right hand side is replaced by i.e. where the collision frequency is constant [18,19]. There are also results concering solutions close to a global equilibrium, which hold also for density and temperature dependent collision frequencies ( [25,26]). There is a rather large litterature concerning various aspects of the BGK-equation dealing, for example, with methods for numerical treatment of rarefied gases (some recent examples are [12,2,24]), their fluid dynamical limits (see for example [20,9,8]), or models accounting for polyatomic gases or mixtures of different gases (for example in [1,11,4]), to give a few examples. A paper attempting to find a well-motivated appoximation of the collision frequenecy 1/σ can be found in [21].
The BGK is a fenomenological equation in the sense that it is derived explicitly to satisfy certain physical properties of a dilute gas, but to our knowledge there are no published works attempting to justify the equations directly from the dynamics of an Nparticle system. This is in contrast with the Boltzmann equation, for which there is now a rigorous derivation starting from the Liouville equation for hard sphere dynamcis, or for short range potentials, of an N -particle system [10,16].
This paper is a step towards providing a justification of the BGK equation. The BGK equation without electric field can be interpreted as a model of a large system of particles, where each particle moves independently with their own (constant) velocity, but where the velocity jumps at exponentially distributed intervals. The jump rate is proportional to the local density of the gas, and after the jump the particle velocity is a normally distributed random variable, independently of the initial velocity, but with mean and variance determined by the local temperature and mean velocity of the gas. One can make a similar interpretation of the Boltzmann equation for hard spheres, but with two important differences. First, in (1), the collision rate only depends on the local density, and not on the velocity of the particles. In fact, equation (1) corresponds to a system of so-called Maxwellian molecueles and not to hard spheres. The second, and more important, difference lies in that the distribution of velocities of particles after a jump. The Boltzmann equation for Maxwellian molecules, which in similar notation is represents a process in which the velocity particle after the jump is given by the outcome of a random collision with a second particle drawn from the distribution with density v → f (x, v, t)/ 3 f (x, w, t) dw. The solutions to (5) converge to a Maxwellian distribution when t → ∞, or equivalently, when the average number of velocity jumps that one particle has made, goes to infinity. In the BGK model, the velocity of a particle has a normal distribution after only one jump, and a particle system that converges to a solution of the BGK model must achieve that in the limit of infinitely many particles.
The particle system that we propose consists of particles that can have two states, active and passive, where only the active particles participate in collisions with other particles. The BGK equation will describe the evolution of passive particles in the limit of infinitely many particles. One may think of the active particles being ions, that interact at a high rate with each other, the passive ones being neutrals that do not interact. On the other hand, a neutral particle and an ion may encounter and interchange state by the transfer of an electron, so that the collection which is similar to allowing the velocity of a neutral particle jump to a random velocity given by the distribution of the acitve particles. And if the collision rate for active particles is very high, then the active particles will have time to come close to an equilibrium distribution before the next exchange with the passive particles takes place.
At a formal level, one may actually pursue these ideas to derive a BGK equation of the form (1), or a hard sphere version of the same, but to make a completely rigorous derivation along the lines of for example [10] seems to be difficult. Long before a rigorous results on the validity of the Boltzmann equation had been obtained for a real particle system, Mark Kac [15] proposed a Markov jump process for the velocities of an N particle distribution, and proved that in the limit as N → ∞, the velocity distribution of one particle converges to the solution of a Boltzmann-like equation for a spatially homogenous gas of Maxwellian molecules with one-dimensional velocities.
In this paper we will construct a Kac-type model of a system of N passive and M active particles, and a jump process involving collisions between active particles and the switch between active and passive state, as described above. We then prove that the one-particle distribution for passive particles converges to a BGK equation of the form where is the standard normal distribution in one dimension. This limit can be obtained in a scaling where M /N → 0 when N → ∞, that is, when the fraction of active particles vanishes in the limit.
The paper is based on the results in the doctoral thesis of the first author [17]. A very similar model, with two different kinds of particles, has been presented by Bonetto et. al. in [6], and also in [5]. The authors are in general interested in kinetic models coupled with a thermostat, and in the cited papers the larger set of particles (N in our paper) is considered as a thermostat acting on the smaller set of particles, and they prove that indeed, when N → ∞ the large N -particle system is a good approximation of a Gaussian thermostat. Related results can also be found in [22].
The paper is organised as follows. In Section 2, we discuss Markov jump processes that give BGK-like equations in the limit of infinitely many particles, and present in full detail our final model. In Section 3 we introduce some further notation, and present the initial steps of the proof. An important step of the proof is to show that the energy partition between the passive and active particles in the limit is such that the mean energy for the active particles is one. This is proven in Section 4 by performing very explicit calculations of certain moments of the solutions. The proof is then concluded in Section 5, which also contains some final remarks.

THE PARTICLE SYSTEM, AND ITS LIMITING KINETIC EQUATION
We consider a particle system consisting of N + M particles, where N is number of passive particles represented by V = (v 1 , . . . , v N ) ∈ N , and M is number of active particles represented by W = (w 1 , . . . , w M ) ∈ M . One active particle is assumed to have the same mass as one passive particle, and here that mass is set to 1. The total kinetic energy, which thefore is 1 2 v 2 1 + · · · + v 2 N + w 2 1 + · · · + w 2 M , is assumed to be conserved, and therfore the state space of the particle system is S N +M −1 ( N + M ), the N + M − 1dimensional sphere of radius N + M . Throughout the paper we also assume that N > M .
The dynamics of the system consist of two independent jump processes. The first one is the so-called Kac walk, which mimics the pairwise collisions of a rarefied gas. This process involves only the active particles. The second process involves a pair consisting of one active and one passive particle, and leads to an exchange of state: the passive particle becomes active, while retaining its velocity, and the active particle becomes passive. In this way there is an exchange of energy between the two sets of particles.
The Kac walk on the set of active particles is defined as follows: • the jumps occur exponentially distributed with rate λ 2 N M .
• in a jump, a pair (w i , w j ) is is chosen uniformly among the active particles, and and θ ∈ [−π, π[ is drawn from an even distribution. Then The jump rate for exchange between active and passive particles is chosen as λ 1 N , which means that the jump rate for a given indexed passive particle v j is λ 1 , independently of N , and that the rate at which active particles become passive is λ 1 N M −1 per active particle. Therefore the subsystem of active particles on average experience λ 2 λ 1 M jumps of Kac-type between two exchange jumps. Without loss of generality we set is λ 1 = 1, and denote the parameter λ 2 simply as λ in what follows.
The intuitive picture is this: consider a time interval [t 1 , t 2 [, where the end points are given by two consecutive exchange events. In this interval the vector V = (v 1 , ..., v N ) is unchanged, and the vector W = (w 1 , ..., w M ) will make on the order of λM N (t 2 − t 1 ) ∼ λM steps in the Kac walk. The energy of the set of active particles is conserved by this process, and hence |W | 2 = w 2 1 + . . . + w 2 M is constant. If λ is very large, the Kac walk will drive the distribution of W to an almost uniform distribution on the sphere defined by |W (t 1 )| 2 . At t 2 a new exchange event takes place, when a randomly chosen active particle becomes passive, and hence the set of passive particles will gain a particle drawn from a distribution which is the marginal of the uniform distribution of an M − 1-dimensional sphere. But this marginal distribution is close to a Gaussian when M is large, and therefore, looking only at the distribution of passive particles, this will loose particles at exponential rate λ 1 = 1 per passive particle, and gain particles drawn from a Gaussian distribution with the same rate; this is the BGK-process for a spatially homogenous gas.
All of this can be quantified, but some notation is needed in order to formulate a theorem. First of all we define the master equation, or forward Kolmogorov equation, corresponding to the jump process. Let F N M (V, W, t) be the probability density with respect to the induced measure σ 1 on S N +M −1 ( N + M ) for the velocities of the particles at time t. The time evolution of F N M is given by the equation: The operator L N M λ defined in (8) is the generator of the original Kac master equation acting on the W -variables, but with a factor λN in front, to give the jump rate as described above. The operator U N M defined in (9), with is the generator of the exchange process, when a passiv and an active particle exchange their state. An essential assumption here, just like in Kac's orginal work, is that F N M is symmetric with respect to permutations of the coordinates of V and of the coordinates of W . This is to say that all passive particles are identical, and identically distributed, and that the same holds for the active particles. Hence any choice of n passive particles is equivalent to choosing the first n. The following notation will be useful: and W m = (w m+1 , ..., w m ) . 1 The symbol σ is used throughout the paper to denote the measure on the sphere S n−1 (r) induced from the Euclidian measure in n , and therefore it is defined only in combination with the domain of integration.
where g(V n ) and h(W m ) are any bounded continuous functions on n , m , respectively, and we assume that f The objective of this paper is to prove that when N , M → ∞ the density of one passive particle, f (10) N M (v, t) converges to a function f (v, t) that satisfies the spatially homogeneous BGK-equation (5). This is formulated in the following theorem: where and where f (v 1 , t) solves the homogeneous BGK equation, So, at least weakly, the one-particle marginal of the N + M dimensional particle system converges to the solution of a BGK equation, as announced in the introduction. The theorem is stated to hold uniformly for t ≥ t 0 > 0, but to achiefe convergence uniformly for all t > 0 one must make stronger assumptions on initial data. If the initial data are chaotic, i.e. meaning that the many-particle marginals are close to products of functions of the coordinates, the L 2 -norm in the theorem grows exponentially in N + M , and choosing C N M ∼ c N +M in equation (13) gives natural class of inital data for which the theorem holds.

INITIAL STEPS OF THE PROOF
To prepare for the proof of Theorem 2.2, we first present a few well-known formulae concerning spheres. A first observation is that although the V -and W -variables represent particles in different states, they behave exactly as variables for integration over the sphere Ω 0,0 = {|V | 2 + |W | 2 = N + M }, and therefore, for any function G(V, W ), we have The area of an n − 1 dimensional sphere of radius r, For any function f defined on S n−1 (r) one may write Therefore the marginals defined in Definition 2.1 may be written explicitly as where we introduce the notation We also define the average with respect to the W -variables as follows: Finally we compute the marginal of the first coordinate of a point chosen uniformly on an M -dimensional sphere of radius M . The uniform density is given by the constant function S M −1 −1 M −(M −1)/2 , and hence the marginal of the first coordinate is The first step in our proof of Theorem 2.2 is to integrate over the W -variables in Equation (7), to find an evolution equation for the V -marginal of F N ,M . It is The integral of L N M λ F N M (V, W, t) vanishes because the generator of the Kac walk conserves mass, and we also use the symmetry with respect to permutations of the Wcoordinates to replace the sum over k in (9) with M terms, all involving w 1 . Adding and subtractingF N M (V j w 1 ) we obtain We will show that by a suitable choice of λ, which is hidden here because it only affects the Kac operator L N ,M ,λ , the term I 2 (V ) vanishes in the limit, and hence that the N M (V ) essentially is governed by I 1 (V ), which in turn will reproduce the righthand side of the BGK-equation in the limit when N , M → ∞.
We denote the N terms of I 1 (V ) as I 1, j (V ), so that I 1 (V ) = N j=1 I 1, j (V ). For j = 1, and for any function g ∈ C( ), by the argument in equation (17), and therefore it is sufficient to consider the first term, I 1,1 (V ). Using Equation 20 and Definition 3.1, the integral in this term is where we have made the change of variables W → τ(V ) W , with the average energy per active particle. The integral is then the marginal distribution of the uniform density over an M − 1-dimensional sphere, as in Equation (23), so that finally the I 1,1 (V ) becomes This integral can now be written as a sum of three terms as follows: Consider the first of these terms, a(V, t). For arbitrary g(v) ∈ C( ), which converges to the right-hand side of equation (14) in Theorem 2.2. Therefore the proof can be concluded by proving that the other terms vanish.
The second term, an integral of b(V, t), converges to zero, because and we know that the M (w) → (w) pointwise, when M → ∞. Of the three terms a(V, t), b(V, t), and c(V, t), the last one is the most difficult to analyse. This is the subject of Section 4, where it is proven that on the domain of integration, we have τ(V ) → 1 when N , M → ∞ under the constraints given in Theorem 2.2.
The proof of Theorem 2.2 can be concluded with these three estimates, together with a proof that g(v 1 )I 2 (V ) d v → 0 when N , M → ∞. The remainig part of this section is devoted to proving that I 2 (V ) converges to zero with suitable choices of N , M , and λ. The result is largely due to Lemma 3.2 below, which in turn follows from a result on the spectral gap for the generator L N M λ of the Kac walk. Without the operator U N M , the generator for the exchange between passive and active particles, equation (7) becomes which is the original Kac master equation in the M variables (w 1 , ..., w M ), with the Vvariables appearing only as parameters. Kac conjectured that the spectral gap ∆ M of − M = −(λN ) −1 L N M λ is bounded away from 0, uniformly in the number of particles M . The parameter V is of course not present i Kac's work, but it doesn't have an influence on the spectral gap, because this gap does not depend on the total energy of the system, and λN only serves to increase the jump rate. Kac's conjecture was first proved by Janvresse, [14], and an exact formula for the gap of − M was later obtained by Carlen, Carvalho and Loss in [7]: and the corresponding eigenfunction is where r is the radius of the M − 1-dimensional sphere. In fact, this eigenvalue and eigenfunction was computed also in [13], but without a proof that this is also determines the spectral gap.
It follows that in an interval t 1 < t < t 2 defined by two consecutive exchange events, With λ very large, the term U N M F can be considered to be a small perturbation, which is expressed in the following lemma: where the norm is in L 2 Ω N ,0 , dσ(W ) .
Proof. By the Duhamel formula we can write The formula for the spectral gap for the Kac model yields A simple computation concludes the proof.
The desired estimate of I 2 (V ) is a direct consequence of Lemma 3.2: W, t) be the solution of equation (7) andF N M (V, t) be given by Definition 3.1 (equation (22)). Then, for every bounded function g : → and all t ≥ 0 (37) Proof. Multiplying equation (7) by F N M (V, W, t) and integrating over Ω 0,0 with respect to σ(V, W ), and using that L M N λ is a non-positive operator, it follows that According to the definition of U N M (equation (9)), the righthand side of this expression is a sum of terms of the form , t) , which due to the inequality 2a b ≤ a 2 +b 2 is smaller than or equal to 1 2 After integration all these terms give a non-positive contribution, and hence the L 2norm of F is non increasing. Therefore the righthand side of the inequality (38) is non-positive, and we have (39) and F N M (V, W, t) ∈ L 2 (Ω 0,0 ) for all t ≥ 0. The function I 2 (V ) is defined as the second sum in the right hand side of equation (25). When multiplying that expression with a function g(v 1 ) (a function depending only on one variable), only the term with j = 1 remains, and therefore Using the Cauchy-Schwartz inequality, and Lemma 3.2 we get A calculation using (9) shows that (40) and so, collecting all the inequalities we finally obtain the inequality (37), which concludes the proof.

EVOLUTION OF MOMENTS AND ENERGY PARTITION
The average energy per active particle is one of them being We will prove that for all t > 0, ψ(t) → 0 when N , M , λ → ∞ in a suitable way.
The starting point is equation (34), where, simplifying notation, F t = F N M (V, W, t), and the operators L = L N ,M ,λ , and U = U N ,M are given in equation (8) and (9). The two operators U and L are self adjoint, but they don't commute. Multiplying the terms in (43) with H(V ) and integrating gives because L acts only in the W -variables, and constants are left invariant by L.
The calculations for moments of the form H(W ), are much simplified when H(W ) is an eigenfunction of the operator L with eigenvalue λ H . In this case In all the integrals we need to compute respectivley.
The mean energy per active particle in a given configuration (V, W ) is and in addition to τ(V ) we introduce the notation as well as the average of m 4 (W ) over the sphere Ω N ,0 . This average can be expressed in terms of the radius of Ω N ,0 , |W |, which in turn is a function of |V |. We have (52) Using this notation, we define the following moments: The expression m 4 (W ) − m 4 (V ) is the eigenfunction of the operator L corresponding to the eigenvalue λN ∆ M . We now obtain expressions for these moments using the equations (44) and (45) 2 .
and therefore, with ν N ,M = 1 + N /M , which means that the mean energy per active particle converges exponentially to 1 as Similarly, to find an expression for ψ(t), we take and obtain the expression 2 Some of the calculations are rather messy, and we have used a computer algebra system for checking these calulations as well as for analysing the linear system for the moments. The code showing all operations are available upon request from the corresponding authors.

Using
After some manipulations, one then gets an equation for ζ(t) Differentiating these expressions, we find a linear system of differential equations for The initial values depend on the moments of the initial density F N M (V, W, 0), and are bounded by These bounds are the worst possible, and obviously are not bounded uniformly in N and M . As we shall see, this is not very critical for η(t), ψ(t), and ζ(t), because after an initial interval of length ∼ M /N , the transient part of the solution will be small for these variables. However, this is not the case for ξ(t), and to conclude our proof of convergence to the BGK-equation, we shall have to assume the initial data F N M (V, W, 0) are such that ξ(0)/M → 0 when M → ∞.
We continue with an asymptotic analysis of the system of equations (65). The matrix A and vectors b 1 and b 1 are expressions involving N and M . We have A = A 0 +A r , where We use Sylvester's formula to compute the exponential exp(tA): where Therefore and this expression can be evaluated at least asymptotically as N /M , M , and λ increase to infinity. For the purpose of this paper we need to prove that for any t > 0 (and uniformly for t ≥ t 0 > 0), ψ(t) → 0 when N /M , M , and λ → ∞ as stated below, but all components of Ψ(t) are needed to obtain a closed system.
The components of the matrices are all rational expressions of N /M , M , and λ∆ M , therefore all terms involving a e ℓ j t , j = 0, 1, 2 vanish when M , N /M and λ∆ M increase to infinity and for t > 0, uniformly for t ≥ t 0 > 0. It is therefore enough to study the terms that are constant in t or with an exponential factor e ℓ 3 t (we recall that ℓ 3 ∼ −1 in the limit of interest).
The matrices A j , j = 1, 2, 3 are, asymptotically, bounded by a constant multiplied by respectively, and imposing that N /M 2 → 0 when N → ∞ the first two terms of (76) vanish in the limit. The condition that N /M 2 → ∞ may be relaxed to N /M 4 → 0 if instead we require ψ(0)/M 2 → 0 when M → 0. And we do need to impose that Summarizing these estimates we obtain the following result: Fix t 0 > 0. Then for t ≥ t 0 there is a constant depending on t 0 such that when N , M , and λ go to infinity.
Remark 4.2. The lemma states that the mean energy per active particle converges to one when the number of particles increase as stated in the lemma. As presented here this is only certain for strictly positive times, but with further assumptions on the initial data F N M (V, W, 0) the same result could be achieved uniformly in time.
We are ready to prove the main result of this section, which says that when N , M , and λ go to infinity as in the previous lemma, the first marginal of from equation (27) converges to zero.
and this converges to zero when N , M , and λ increase to infinity according to Lemma 4.1.
Proof. Multiplying c(V ) with g(v 1 ) and integrating gives This integral may be estimated by noting that (23)), and similarly when M is large enough. In equation (82) we may then estimate g(v 1 ) with g ∞ , and then carry out the integral over u to get The constant C here is obtained by integrating the expression (83). Then Lemma 4.1 provides the needed bounds for ψ(t)

PROOF AND CONCLUSIONS
The goal of this section is to prove Theorem 2. and hence that, for any bounded, continuous function g(v 1 ), These terms have been analysed above in this paper, and it only remains to put the pieces together. The first term, the integral of a(V ), converges to the righthand side of equation (14) ( see eq. (28)), The second term, the integral of b(V ), converges to zero because M (w) converges pointwise to the Maxwellian (see eq. (29)). That the third term converges to zero is exactly the content of Lemma 4.3, and finally Lemma 3.3 states that converges to zero if F N M (·, ·, 0) L 2 (Ω 0,0 ) /λ → 0. By hypothesis (see equation (13)), F N M (·, ·, 0) L 2 (Ω 0,0 ) ≤ C N M for a family of constansts C N M and we may the choose λ = λ N M accordingly. . Now let f (v, t) be the solution of Then and it follows that converges to zero under the assumptions of the theorem.
This concludes the proof of Theorem 2.2. We end by two remarks. First, the theorem is stated to hold for t > t 0 with t 0 > 0 arbitrary. In fact, as stated the theorem is true for any t > 0, but unless the initial data is well prepared the rate of convergence is not necessarily uniform as t → 0.
The second remark concerns the propagation of chaos, i.e. the property that if any m-particle marginal of the initial data factorizes as a product of m one-particle densities when N , M go to infinity, then same holds for the m-particle marginals of the solutions F N M (V, W ). In our case only marginals for the components of V are interesting, and then the result follows rather easily because the jumps never involve more than one component of V at the same time, and because of the structure of the limiting equation.