A Consistent BGK Model with Velocity-Dependent Collision Frequency for Gas Mixtures

We derive a multi-species BGK model with velocity-dependent collision frequency for a non-reactive, multi-component gas mixture. The model is derived by minimizing a weighted entropy under the constraint that the number of particles of each species, total momentum, and total energy are conserved. We prove that this minimization problem admits a unique solution for very general collision frequencies. Moreover, we prove that the model satisfies an H-Theorem and characterize the form of equilibrium.


Introduction
In this paper, we present a BGK-type model for gas mixtures that, in the case of two species, takes the form along with appropriate boundary and initial conditions. Here f 1 = f 1 (x, v, t) and f 2 = f 2 (x, v, t) are the number densities of species of mass m 1 and m 2 , respectively, with respect to the phase space measure dxdv; x ∈ R 3 is the position coordinate of phase space; v ∈ R 3 is the velocity coordinate; and t ≥ 0 is time. The relaxation operator on the right hand side of (1) involves target functions of the form which depend on parameters λ k j = (λ , and (non-negative) collision frequencies ν k j . These parameters depend implicitly on f 1 and f 2 , and once specified, determine the BGK operator.
The purpose of the relaxation operator in (1) is to provide an approximation of the multi-species Boltzmann collision operator that is more computationally tractable, but still maintains important structural properties. In the single-species case, the original BGK model [2] serves this purpose. In particular, it has the same collision invariants as the Boltzmann operator (which lead to conservation of number, momentum, and energy) and it satisfies an H-Theorem. In the multi-species case, these requirements are not as straight-forward to satisfy, but it can be done. There are many BGK models for gas mixtures proposed in the literature [1,5,10,12,[14][15][16]23,28], many of which satisfy these basic requirements and, in addition, are able to match some prescribed relaxation rates and/or transport coefficients that come from more complicated physics models or from experiment. Many of these approaches have been extended to accommodate ellipsoid statistical (ES-BGK) models, polyatomic molecules, chemical reactions or quantum gases; see for example [3,4,13,[24][25][26][27]31].
A common feature of all the models mentioned above is that they only allow for collision frequencies which are independent of the microscopic velocity v of the particles [30]. However, the collision frequencies in principle should depend on the microscopic velocity, which is typically neglected for the reason of simplicity. In the case of neutral gases, velocity independent collision frequency leads to transport properties in the fluid regime that are inconsistent with the full kinetic collision operator, e.g., the Prandtl number. Models such as the ES-BGK model and the Shakov model make changes to the target Maxwellian to provide extra degrees of freedom to the system and correct this defect, but still retain the constant collision frequency assumption. Some attempts have been proposed to re-introduce velocity dependence in the case of variable hard spheres interactions for neutral gases [22], for which velocity-dependent collision frequencies are monotonically increasing and are well-defined.
In the case of Coulomb interactions, i.e. plasmas, the Boltzmann cross section is even more strongly dependent on the relative velocity of the particles, as particles with small relative velocity, or grazing collisions, are the dominant event. As a result, the cross section is nearsingular in both the relative velocity and scattering angle variables, and approximating it with a single, constant value is likely to miss the rich underlying features of the grazing collision physics. Indeed, in a widely used model in hydrodynamic modeling of plasmas, based on a non-conservative but velocity-dependent BGK model [21], the exponent of the velocity component directly appears in the formulas for electrical and thermal conductivity. While we derive a conservative velocity dependent BGK model in this paper, we follow the spirit of [21] and define our velocity-dependent collision frequency in terms of the momentum transfer cross section, which results in a collision frequency that is decreasing in the limit of large relative velocities [20].
The velocity-dependent BGK model of this paper provides two advantages over more accurate Boltzmann and Fokker-Planck-Landau collision models from a computational perspective. Though computation of this model is not the focus of this paper, we mention these motivating concerns here for completeness; numerical discretization will be the subject of future work. First, like other BGK models it is amenable to implicit time discretization, which allows for time steps that exceed the frequently stiff collision frequency. While some proposed methods allow for time steps larger than the collision frequency in the Boltzmann or Landau-Fokker-Planck models [17], these methods are only weakly asymptotic preserving and it can be unclear without a priori knowledge of the size of the deviation from local Maxwellians. Furthermore, multispecies extension of these methods [18] penalize each species with a single relaxation operator rather than penalizing individual reaction pairs, which can lead to inaccuracies in cases where some collisional combinations are more important than others.
Another concern is the computational cost of a collision operator evaluation. The main cost of the BGK model presented here lies in the evaluation of the parameters in the target Maxwellians. When the collision frequencies are constant with respect to velocity, these parameters can be computed explicitly in terms of the given moments. In the velocitydependent case, an optimization problem (similar to (8) and (15) presented below) must be solved numerically via a dual formulation. The main cost here lies in a Newton solvespecifically, the evaluation of the integrals that appear in the gradient and Hessian of the dual objective functions. In practice, the cost of the optimization will depend on the details of the implementation [34][35][36][37]. However, assuming that the velocity grid that is used to discretize the BGK equation in velocity space is also used to evaluate the integrals, then the method will scale like O(N 3 ). For comparison, currently the fastest algorithms for the evaluating the Boltzmann collision operator are spectral methods. In three velocity dimensions, the complexity of evaluating specialized kernels is O(M N 3 log N ) [32], while for general kernels, the complexity is O(M N 4 log N ) [33]. Here N is the number of points in each dimension of the velocity space and M is the number of quadrature points needed to accurately approximate integrals over the unit sphere S 2 in R 3 . In practice, the size of M is problem dependent; and while M N 2 is typical, it is not unusual to have M ≥ N [33]. Thus, while more expensive than the constant frequency case, the BGK model with velocity-dependent frequencies is still expected to be less expensive than evaluating the Boltzmann collision operator.
In this paper, we derive a model of the form (1) that allows for velocity-dependent collision frequencies. Our derivation includes as a by-product the single-species BGK model with velocity-dependent collision frequency that was proposed in [29]. We identify target functions that are consistent with the conservation laws for (1) and satisfy an entropy minimization principle. In particular, intra-species collisions (between the same species) should preserve mass, momentum, and energy within a species; that is, Meanwhile inter-species collisions (between different species) should preserve the mass of each species, but only the combined momentum and energy of both; that is, When the collision frequencies are independent of v, the integrals in (3) and (4) can be computed explicitly, thereby providing relationships between the parameters λ k j and the moments of f 1 and f 2 with respect to {1, v, |v| 2 }. In the single-species case, this relationship defines the target function as the Maxwellian associated to f , while in the multi-species case, additional constraints must be imposed. However, when the collision frequencies depend on v, the aforementioned integrals are not always computable in closed form and the relationship between the parameters λ k j and the moments of f 1 and f 2 with respect to {1, v, |v| 2 } cannot be written down analytically.
In spite of the difficulty of relating the target parameters to the moments of the kinetic distributions, the entropy minimization formulation can be still used to establish a unique set of parameters, under the conditions λ 12 1 = λ 21 1 and λ 12 2 = λ 21 2 . We do so by adapting the strategy from [19] to fit the current setting. While a more abstract approach based solely on convex optimization tools can also be used [6], we follow [19] because it provides a more concrete connection to the application at hand. Our proof provides a rigorous justification for the target function used in [29] for the single species case. It also leads to an H-Theorem for the multi-species system (1).
The remainder of the paper is organized as follows. In Sect. 2, we motivate the choice of the target Maxwellians as solutions of minimization problems of the entropy under certain constraints. In Sect. 3, we prove existence and uniqueness of the minimization problems. In Sect. 4, we prove consistency of the model meaning that it satisfies the conservation properties, the H-Theorem and Maxwell distributions with equal mean velocity and temperature in equilibrium. In Sect. 5, we briefly summarize the straightforward extension to the case of N species, still with binary interactions.

The Structure of the Target Functions
In this section, we motivate the form of the target functions in (2). It will be convenient in what follows to define the strictly convex function and the vector-valued function Since h is convex and h (z) = ln(z), it follows that

The One Species Target Maxwellians
We seek a solution of the weighted entropy minimization problem where The choice of the set χ k ensures the conservation properties (3) for intra-species collisions.
The motivation for weighting the usual objective by the collision frequencies in (8) is that the ansatz will take the form (2). Indeed, by standard optimization theory, any critical point satisfies the first-order optimality condition which implies then that In Sect. 3.1, we prove in a rigorous way that there exists a unique function of the form (12) that satisfies these constraints.

Theorem 1
Suppose that there exists λ kk ∈ R × R 3 × R such that the function M kk given in (12) is an element of χ k . Then M kk is the unique minimizer of (8).
Proof According to (7) point-wise in v. Thus, because ν kk ≥ 0, it follows that for all g ∈ χ k , Hence M kk is a minimizer of (8), and uniqueness follows directly from the strict convexity of h.

The Mixture Target Maxwellians
For interactions between species, we seek a solution of the weighted entropy minimization problem min g 1 ,g 2 ∈χ 12 where Here, χ 12 is chosen such that the constraints (3) for inter-species collisions are satisfied. Similar to the case of intra-species collisions, we consider the Lagrange functional L : Any critical point (M 12 , M 21 , λ 1 0 , λ 2 0 , λ 1 , λ 2 ) of L satisfies the first-order optimality conditions where λ 12 = (λ 1 0 , λ 1 , λ 2 ) and λ 21 = (λ 2 0 , λ 1 , λ 2 ). Therefore Since we only require conservation of the combined momentum and kinetic energy, there is only one Lagrange multiplier for the momentum constraint and one Lagrange multiplier for the energy constraint. Therefore, λ 12 1 = λ 21 1 and λ 12 2 = λ 21 2 in (2). When the collision frequency is constant, this restriction is the same as the one used in [15], but more restrictive than the model in [23].
In the next section, we prove the existence of functions of the form (2) that satisfy the constraints in (3) and (4). As in the single species case, it follows that these functions are unique minimizer of the corresponding minimization problem.
Proof According to (7) point-wise in v, for any measurable function g and k, j ∈ {1, 2}. Therefore, since ν k j ≥ 0, it follows that for any measureable functions g 1 and g 2 , Since λ 12 1 = λ 21 1 and λ 12 2 = λ 21 2 , If (g 1 , g 2 ) and (M 12 , M 21 ) are elements of χ 12 , then the constraints in (16) imply that each of the terms above is zero. In such cases, (23) reduces which shows that (M 12 , M 21 ) solves (15). Since the collision frequencies ν 12 and ν 21 are non-negative and h is strictly convex, it follows that this solution is unique.

Existence and Uniqueness of the Target Maxwellians
In this section, we prove the existence of the multipliers λ 11 , λ 22 , λ 12 and λ 21 such that the single-species targets M 11 and M 22 satisfy (3) and the mixture targets M 12 and M 21 satisfy (4). We follow closey the strategy laid out in [19], although some variations will be needed to account for the velocity-dependent collision frequencies and the mixture targets. Throughout the paper, we denote a distribution function of exponential form by and let For any g ∈ D k j the moment map μ k j is given by We make the following assumptions about the collision frequencies.
Assumption 1 Each frequency ν k j is strictly positive and defined such that is independent of k and j.
Roughly speaking, these assumptions are used to ensure integrability properties that are satisfied when the collision frequencies are independent of the velocity. They are used in the technical details of the proofs below, but are in practice satisfied by many realistic frequency models.

Target Functions for Intra-species Collisions
We start the intra-species case; that is, for k ∈ {1, 2}, we show the existence of multiplier λ kk such that M kk satisfies (3). The basic idea is to show that the dual function is differentiable and attains its minimum on Λ for any ρ ∈ μ kk (D kk ). Then the necessary condition for an extremum in Λ yields which gives ρ = μ kk (exp k λ kk ).
For any nonzero δ ∈ R 5 where A Taylor series expansion shows that Because Λ is open, for |δ| sufficiently small, exp λ/2+δ (v) and exp λ/2−δ (v) are elements of D kk , in which case g δ is integrable. Moreover, exp λ/2 is bounded. Hence f δ is bounded above by an integrable function and the dominated convergence theorem gives The existence of the Hessian can be proven in an analogous way.

attains its unique minimum in the open interval
takes the value +∞ if the boundary ∂Λ is not met in the direction ξ.
Proof The fact that z is strictly convex and differentiable with respect to λ implies that z ξ is strictly convex and differentiable with respect to s. Hence it attains a unique minimum on the closure of I (ξ, λ). We now show that z ξ cannot attain its minimum on the boundary of I (ξ, λ). Suppose first that s b (ξ, λ) < ∞. According to Assumption 1, λ + s b (ξ, λ)ξ / ∈ Λ. Hence by Fatou's Lemma, which implies that lim s→s b (ξ,λ) z ξ (s) = +∞. Suppose now that s b (ξ, λ) = ∞. There are two cases: , there exists g ∈ D kk such that ρ = μ kk (g). By definition, g is not identically zero and by Assumption 1 ν kk > 0. Thus the set has positive measure. Hence due to exponential growth in s.
Corollary 1 Given any f k ∈ D kk , there exists a unique multiplier λ kk such that M kk given by (2) solves (8).

Lemma 3 The functionz defined in (50) is strictly convex and twice Fréchet differentiable onΛ.
Proof Differentiability of thez can be deduced as in the intra-species case by simply following the arguments of Lemma 1. We skip these details. Convexity also follows in a similar way. Letφ(λ) = μ 12 0 (exp 1 λ 1 ) + μ 21 0 (exp 2 λ 2 ), then convexity of the exponential function implies that for any θ ∈ (0, 1), λ ∈Λ, and β ∈Λ, Thusφ is strictly convex, as isz, since the two functions differ only by a linear term.
Proof We follow the arguments of the proof of Lemma 2. The fact thatz is strictly convex and differentiable with respect to λ implies thatz ξ is strictly convex and differentiable with respect to s. Hencez ξ attains a unique minimum on the closure ofĪ (ξ, λ). We therefore need only show thatz ξ cannot attain its minimum on the boundary ofĪ (ξ, λ).
The proof of this theorem is analogous to the proof of Theorem 3 in the intra-species case.

Consistency of the Model
The conditions (3) and (4) lead to standard conservation laws and an entropy dissipation statement. We recall a few definitions: The mass density, momentum, and energy of an integrable distribution g = g(v) of particles with mass m are given by the moments respectively. The associated mean velocity and temperature are given by

Conservation Properties
An immediate consequence of (3) and (4) is the following.
Theorem 5 (Conservation of the number of each species, total momentum and total energy) The space-homogeneous form of (1) satisfies

Entropy Dissipation and the Structure of Equilibria
Define the total entropy density and the dissipation density = ν 11 ln g 1 (M 11 − g 1 )dv + ν 12 ln g 1 (M 12 − g 1 )dv (73) Theorem 6 Assume g 1 , g 2 > 0. Then S(g 1 , g 2 ) ≥ 0 with equality if and only if g 1 and g 2 are two Maxwellian distributions with equal mean velocity and temperature.
Proof In [29], it is shown that S kk (g) ≥ 0 with equality if and only if g is a Maxwellian. Thus it remains to show a similar result for the combined quantity S 12 (g 1 , g 2 ) + S 21 (g 1 , g 2 ).
We begin with the following claim: Indeed an explicit calculation gives which when substituted into (75) gives due to the constraints (4). From (75), it follows that with equality if and only if g 1 = M 12 and g 2 = M 21 . Moreover, a direct calculation shows that the functions M 12 and M 21 have the same mean velocity and temperature: Corollary 3 (Entropy inequality for mixtures) Assume that f 1 , f 2 > 0 are a solution to (1) where the target Maxwellians have the shape (2), then we have the following entropy inequality with equality if and only if f 1 and f 2 are two Maxwellian distributions with equal mean velocity and temperature.
Proof A direct calculation with (1) gives The result then follows immediately from the previous theorem.

The N-Species Case
The two-species case can be extended to a system of N -species that undergo binary collisions. We consider the N -species kinetic equation, The quantity ν ii is the collision frequency of particles of species i with itself whereas ν i j is the collision frequency of particles of species i with species j, with i, j = 1, . . . , N , i = j.
We only have terms of this form and not terms containing indices of more than two species because we consider only binary interactions. For fixed i, j ∈ {1, . . . , N } the target Maxwellians M ii , M j j , M i j and M ji are given by (2). The single species target Maxwellians M ii and M j j will be determined such that they satisfy (3). The functions M i j and M ji will be determined such that we obtain conservation of mass of each species and conservation of total momentum and total energy in interactions between these two species, i.e., as an obvious generalization of (4). All the proofs concerning existence and uniqueness of the target Maxwellians and the H-Theorem can be proven exactly in the same way as for two species. For the total entropy H ( f 1 , . . . , f N ) = (h( f 1 ) + · · · + h( f N ))dv we obtain ∂ t (H ( f 1 , . . . , f N )) + ∇ x · v(h( f 1 ) + · · · + h( f N ))dv ≤ 0.

Conclusion
We have presented a multi-species BGK model in which the collision frequencies depend on the microscopic velocity. The model is formally derived based on an entropy minimization principle, which implies that the target functions take the form of Maxwellians. However, contrary to classical BGK models with velocity-independent frequencies, the relationship between the Maxwellian parameters and the moments of the distribution function is not analytic. Thus some effort is required to establish rigorously the existence of parameters which satisfy first-order optimality conditions. We also show that the derived model satisfies an H-Theorem and that it can be extended to the case of arbitrarily many species undergoing binary collisions.
In future work, we will develop numerical tools for discretizing the model developed here, including the numerical solution of the defining optimization problem. A numerical code will enable computational explorations about how to choose the collision frequencies and what benefit is providing by their flexibility. Also, because the motivation for the model is the simulation of multi-species plasmas, we will extend it for use in such contexts by adding self-consistent fields.