The Boltzmann equation for plane Couette flow

In the paper, we study the plane Couette flow of a rarefied gas between two parallel infinite plates at $y=\pm L$ moving relative to each other with opposite velocities $(\pm \alpha L,0,0)$ along the $x$-direction. Assuming that the stationary state takes the specific form of $F(y,v_x-\alpha y,v_y,v_z)$ with the $x$-component of the molecular velocity sheared linearly along the $y$-direction, such steady flow is governed by a boundary value problem on a steady nonlinear Boltzmann equation driven by an external shear force under the homogeneous non-moving diffuse reflection boundary condition. In case of the Maxwell molecule collisions, we establish the existence of spatially inhomogeneous non-equilibrium stationary solutions to the steady problem for any small enough shear rate $\alpha>0$ via an elaborate perturbation approach using Caflisch's decomposition together with Guo's $L^\infty\cap L^2$ theory. The result indicates the polynomial tail at large velocities for the stationary distribution. Moreover, the large time asymptotic stability of the stationary solution with an exponential convergence is also obtained and as a consequence the nonnegativity of the steady profile is justified.

The steady state of a rarefied gas between two parallel plates with the same temperatures and opposite velocities is one of the most fundamental boundary-value problems in kinetic theory, see the books of Kogan [28], Cercignani [11], Garzó-Santos [22], and Sone [33]. In particular, numerical analysis of the plane Couette flow of rarefied gas on the basis of the nonlinear Boltzmann equation has been extensively conducted in the physical literatures, cf. [29][30][31]34,35]. On the other hand, the mathematical study on this problem, even in the case when there is a temperature gap between two plates and a constant external force parallel to the boundaries, has been carried out by Esposito-Lebowitz-Marra [18,19] for the hydrodynamic description of the steady rarefied gas flow via the approximation of the corresponding compressible Navier-Stokes equations with no-slip boundary condition. The result in [19] for hard sphere model was later extended in [13] to the case of hard intermolecular potentials with Grad's angular cutoff as well as the Maxwell molecule case for which only the polynomial decay of the stationary solution for large velocities is obtained compared to the exponential decay for hard sphere model. In addition, closely related to the plane Couette flow, the stationary Boltzmann equation for rarefied gas in a Couette flow setting between two coaxial rotating cylinders was also studied extensively by Arkeryd-Nouri [2,3] in the fluid dynamic regime, see also a recent work [1] for further investigation of ghost effect induced by curvature.
The current study of the plane Couette flow with boundaries is motivated by the previous work [14] by the first two authors for uniform shear flow via the Boltzmann equation without boundaries. We refer readers to [9,12,21,27,36,37] and references therein for more details of the topic on uniform shear flow. In particular, in a recent significant progress [9], Bobylev-Nota-Velázquez studied the self-similar asymptotics of solutions in large time for the Boltzmann equation with a general deformation of small strength and they also showed that the self-similar profile can have the finite polynomial moments of higher order as long as the deformation strength is smaller. In this paper, we will take into account the effect of shear force induced by the relative motion of the boundaries. We hope that the current study can shed some light on the relation between the Couette flow with boundary and the uniform shear flow without boundary. A rigorous justification of the behavior of solutions in the limit L → ∞ is left for future research.
To specify the problem, we consider the rarefied gas between two parallel infinite plates with the same uniform temperature T 0 > 0, one at y = +L is moving with velocity (U + , 0, 0) and U + = αL and the other at y = −L is moving with velocity (U − , 0, 0) and U − = −αL, where α > 0 is a parameter for the shear rate, see Figure 1 below. Moreover, we assume that the gas molecules are of the Maxwellian type and reflected diffusively on the plates y = ±L. Let the non-negative unknown function F = F (y, v) ≥ 0 be the time-independent density distribution function of gas particles with velocity v = (v x , v y , v z ) ∈ R 3 located at position y ∈ (−L, L) along the vertical direction with the slab symmetry in the horizontal (x, z)-plane in space. Then, the motion of such rarefied gas flow can be governed by the steady Boltzmann equation For the Maxwell molecule model, the collision operator operator Q, which is bilinear and acts only on velocity variable, takes the form of where the velocity pairs (v * , v) and (v ′ * , v ′ ) satisfy the relation v ′ denoting the ω-representation according to conservation of momentum and energy in the elastic collision, i.e., v * +v = v ′ * +v ′ and |v * | 2 +|v| 2 = |v ′ * | 2 +|v ′ | 2 , respectively. Throughout the paper, we assume that the collision kernel B 0 (cos θ) with cos θ = (v − v * ) · ω/|v − v * |, depending only on the angle θ between the relative velocity v − v * and ω, satisfies the Grad's angular cutoff assumption 0 ≤ B 0 (cos θ) ≤ C| cos θ|, (1.6) for a generic constant C > 0.
In the paper, for the boundary-value problem (1.1), (1.2) and (1.3) with finite Knudsen number, we look for stationary solutions of the following specific form where the horizontal molecular velocity v x − αy is sheared linearly along the y-direction. After plugging (1.7) into (1.1), (1.2), (1.3) and normalizing L, M and T 0 to be one for brevity, the stationary distribution function F st is determined by the following boundary-value problem with the global Maxwellian µ = (2π) −3/2 e −|v| 2 /2 . This paper aims to establish the existence of solutions to the above boundary value problem (1.8) for any small enough shear rate α > 0, as well as its large time asymptotic stability. To solve (1.8), we will apply the perturbation approach by taking the shear rate as a small parameter. If α = 0, F st = µ is the unique equilibrium solution to the boundary value problem (1.8). However, for α > 0, the external shear force drives the rarefied gas far from the equilibrium. Precisely, we set (1.10) By plugging (1.9) into (1.8) and comparing coefficients of the equation in the order of α, we obtain the equation for G 1 v y ∂ y G 1 + LG 1 = −v x v y √ µ, (1.11) with boundary condition √ µG 1 (±1, v)|v y | dv, (1.12) and the equation for the remainder (1.13) with boundary condition (1.14) Here, the linear and nonlinear collision operator L and Γ are given by and respectively. Properties of these two operators will be presented in Section 2. Note that to solve G 1 , both (1.11) and (1.12) with the restriction Hence, the diffuse reflection boundary condition (1.12) for G 1 can be reduced to the homogeneous inflow boundary condition (1. 16) The first result in the paper for the existence of the Couette flow problem is stated as follows. To the end, we use a velocity weight function w q = w q (v) := (1 + |v| 2 ) q (1.17) with an integer q > 0.
Theorem 1.1. Assume that the Boltzmann collision kernel is of the Maxwell molecule type (1.6). Then, the boundary value problem (1.8) admits a unique steady solution F st = F st (y, v) ≥ 0 of the form (1.9) satisfying (1.10) and the following estimates on G 1 and G R , respectively.
(i) The first-order correction G 1 = G 1 (y, v), uniquely solved by the boundary value problem (1.11) and (1.16), satisfies (1.15) and for any integers m ≥ 0 and q ≥ 0, whereC 1 > 0 is a constant depending only on m and q. (ii) The remainder G R = G R (y, v), uniquely solved by the boundary value problem (1.13) and (1.14), satisfies that there is an integer q 0 > 0 such that for any integer q ≥ q 0 , there is α 0 = α 0 (q) > 0 depending on q such that for any α ∈ (0, α 0 ) and any integer m ≥ 0 G R := √ µG R satisfies that 19) whereC m,q > 0 is a constant depending only on m and q but independent of α.
Some remarks on Theorem 1.1 are given as follows.
Remark 1.1. The steady solution F st to the boundary value problem (1.8) is essentially constructed in the regime where the collision is dominated and the shearing effect is weak. By (1. 18) and (1.19), the steady solution takes the form of with the remainder of the second order decaying in large velocities only polynomially. The order of the polynomial decay can be arbitrarily large as long as the shear rate is sufficiently small. It generally holds α 0 (q) → 0 as q → ∞, and in particular, one may take α 0 (q) = ν 0 8q as shown in the proof. The result is consistent with the one in [14] for uniform shear flow without boundaries in the spatially homogeneous setting.
Remark 1.2. Without using the odd-in-v x property as in (1.15), the existence of G 1 (y, v) to the BVP (1.11) under the diffuse reflection boundary condition (1.12) also can be established by the same approach as for treating the remainder G R . Here, we take this formulation only for brevity of presentation because the proof for the homogeneous inflow boundary is relatively easier than that for the diffuse reflection boundary. Remark 1.3. We notice that it is necessary to deal with the v x -derivative estimates due to the appearance of the shear force term v y ∂ vx F st , in particular, the term v y ∂ vx G 1 becomes a source term in the equation (1.13) for G R . We emphasize that although one can obtain the derivative estimates as in (1.18) and (1.19) in v x , it is impossible to obtain similar estimate on derivative in v y because G 1 (y, v) is discontinuous at v y = 0, see (4.25) for an explicit form of G 1 when the non-local collision term is omitted.
To establish the nonnegativity of the stationary profile F st (y, v), we further study the following initial boundary value problem of the Boltzmann equation with a shear force (1.21) One may expect that the solution of the time-dependent problem (1.21) tends in large time toward that of the steady problem (1.8). For this, the second result is concerned with the large time asymptotic stability of the stationary solution F st which gives the nonnegativity of F st .
Theorem 1.2. Let F st (y, v) be the steady state obtained in Theorem 1.1 corresponding to a shear rate α ∈ (0, α 0 ). There are constants ε 0 > 0, λ 0 > 0 and C > 0, independent of α, such that if initial data F 0 (y, v) ≥ 0 satisfy then the initial boundary value problem (1.21) admits a unique solution F (t, y, v) ≥ 0 satisfying the following decay estimate:  ) is uniform in all α ∈ (0, α 0 ) when the large enough integer q is chosen and hence α 0 = α 0 (q) > 0 is fixed. Thus, the exponential time decay estimate (1.23) also holds uniformly for any α ∈ (0, α 0 ), in particular, C and λ 0 are independent of α. As α → 0, we are able to recover the exponential convergence of the solution F (t, y, v) to the global Maxwellian µ in L ∞ -norm weighted by the polynomial velocity weight w q (v).
In what follows we present key points and strategy in the proof of the main results stated above. As pointed out in a recent nice survey by Esposito-Marra [20], stationary non-equilibrium solutions to the Boltzmann equation, despite their relevance in applications, are much less studied than time-dependent solutions, and no general existence theory is available, due to technical difficulties. Readers may refer to [20] and references therein for a thorough review on this subject. As for the Boltzmann equation on the plane Couette flow, [19] and [13] mentioned before seem to be the only mathematical works on the fluid dynamic approximation solutions in the steady case for small Knudsen number. But it remains unsolved how to justify the large time asymptotics toward the stationary solution for the time-dependent problem in the same setting of the fluid limit. In this paper, motivated by [14], instead of constructing the fluid dynamic approximation solutions, we focus on the existence and dynamical stability of the plane Couette flow with the finite Knudsen number for both the steady and unsteady problems.
First of all, for the original Couette flow problem (1.1), (1.2) and (1.3), we note that a direct perturbation approach by linearization of the boundary condition in α in terms of the techniques in [13,18,19] or [16] can be applied to prove the existence of stationary solutions, because the inhomogeneous data appear only on the tangent (x, z)-plane. The solution thus obtained has the structure around global Maxwellians of the form corresponding to the linearization of the wall Maxwellians at y = ±L On the other hand, in the formulation used in this paper, we rather look for the solution of the specific structure (1.7), and hence the problem can be reduced to solve (1.8) for the Boltzmann equation driven by an external shear force under the homogeneous non-moving diffuse reflection boundary condition. This means that the solution to the Couette flow problem (1.1), (1.2) and (1.3) is established around the local Maxwellian µ(v x − αy, v y , v z ) instead of the global Maxwellian µ such that the kinetic diffusive reflection boundary condition (1.2) is satisfied for the background solution µ(v x − αy, v y , v z ). In addition, as mentioned before, it seems more convenient to use the formulation with shear forces than the original one driven by the relative motion of boundaries in order to understand the asymptotic behavior of solutions in the limit L → ∞, that is, how the Couette flow with boundaries converges to a shear flow without boundary that is closely related to what has been studied in the previous works [14] for uniform shear flow in the spatially homogeneous setting.
We also comment on the boundary value problem (1.11) and (1.12) for determing the first order correction term G 1 (y, v). Notice that the inhomogeneous source term −v x v y √ µ in (1.11) does not satisfy the boundary condition (1.12), so a space-dependent non-trivial solution is induced. If the boundary condition is omitted and only the spatially homogeneous equation is considered, the corresponding solution can be written as The form (1.24) is then consistent with the uniform shear flow in [14]. To solve the boundary value problem (1.11) and (1.12), the same approach as for treating the remainder G R can be applied. However, in order to simplify the proof, we have made use of an additional property (1.15) to reduce the diffusive reflection boundary condition (1.12) to the homogeneous inflow boundary condition (1.16). To treat (1.11) and (1.16), we develop a direct L ∞ -L 2 method without using the stochastic cycles as in [24]. In particular, thanks to the splitting L = ν 0 − K, if the non-local term KG 1 is omitted, the solution to the boundary value problem can be explicitly expressed as Moreover, we use the bootstrap argument as in [15] to treat the following problem with a parameter With the solvability starting from σ = 0, we are able to iteratively solve the above boundary value problem for σ over the intervals [0, σ * ], [σ * , 2σ * ], and so on, where σ * > 0 is small enough such that σ * KG 1 can be regarded as a source term in the L ∞ estimation. Therefore, in the end, the original problem corresponding to σ = 1 can be solved. In this procedure, the uniform L ∞ estimate can be obtained through the interplay with the L 2 estimates in terms of the Guo's technique in [24]. Here, we have omitted the discussions on the mass conservation (1.10) for G 1 . In fact, inspired by [24], an extra damping term ǫG 1 with the vanishing parameter ǫ > 0 has to be used, cf. Section 4 for details.
We now discuss some key points about estimating the remainder G R for the existence of solutions to the boundary value problem (1.13) and (1.14). The direct L ∞ -L 2 approach is no longer available because the linear term 1 2 αv x v y G R can not be controlled in the large velocity regime. Notice that this term arises from the action of the shear force on the exponential weight function √ µ in the perturbation. To overcome it, as in [14], we apply the Caflisch's decomposition and G R,1 and G R,2 satisfy the coupled boundary value problems respectively. Then, in (1.25), the term − 1 2 α √ µv x v y G R,2 can be controlled due to the appearance of √ µ. Here, since the operator norm of K may not be small, the term χ M KG R,1 over the large velocity regime can be viewed as a source in (1.25) for G R,1 , while the complementary term (1 − χ M )µ − 1 2 KG R,1 is taken as a source in (1.26) for G R,2 . A crucial observation inspired by [4] in estimating G R,1 is that the norm of the weighted operator w q χ M K on L ∞ v with the polynomial velocity weight w q = (1 + |v| 2 ) q can be arbitrarily small as long as M and q are chosen sufficiently large, see Lemma 2.4. Notice that Lemma 2.4 holds only for the Maxwell molecule potential as shown in the proof. Compared to the previous work [14] for uniform shear flow, it is more complicated to solve the coupling steady boundary value problems (1.25) and (1.26) in a bounded domain. We now list the main steps in the proof.
• Step 1. We first modify the coupled boundary value problems with two parameters ǫ > 0 being small enough and 0 ≤ σ ≤ 1, see (5.10), and obtain the a priori estimates uniform in ǫ and σ in the L ∞ framework, see Lemma 5.1 and the proof for Proposition 5.1. For the proof of Lemma 5.1, we apply Guo's approach in [24] to the shear flow problem in a slab. In particular, we introduce the mild formulation (5.23) to treat the diffuse boundary condition with the help of Lemma 8.1, and reprove Ukai's trace theorem in Lemma 3.1 for the L 2 estimates. • Step 2. Similar to solving the first order correction term G 1 , we design an explicit procedure to solve the parametrized boundary value problem (5.10) iteratively for σ ∈ [0, 1] from σ = 0 to σ = 1 for any fixed ǫ > 0, see Lemma 5.2. Notice that the problem for σ = 0 is reduced to the one without the nonlocal collision terms under the homogeneous inflow boundary condition so that the characteristic method can be directly applied.
• Step 3. We study the limit ǫ → 0 to obtain the desired solution, see Subsection 5.4 for details. The key point is to obtain the macroscopic estimates in order to bound the L 2 norm of G R,2 in terms of the L ∞ norm of G R,1 . We apply the dual argument developed first in [16]. Note that it is delicate to deduce these estimates to be uniform for any small parameter ǫ > 0.
With the existence of stationary solution F st , the asymptotic stability of the perturbation F = F st + √ µf as (6.1) is considered in the reformulated IBVP as (6.2). Technically, we follow the same strategy as for treating the steady problem. Precisely, we also use the decomposition √ µf = f 1 + √ µf 2 with f 1 , f 2 satisfying the coupled IBVPs (6.4), (6.5) and (6.6), (6.7), respectively. In order to treat initial data with only the polynomial velocity weight, we set f 2 (0, y, v) ≡ 0 and the boundary conditions on f 1 and f 2 both as diffuse reflections which are slightly different from (1.25) and (1.26) in the steady problem. Moreover, in contrast with the steady case, we need to construct suitable temporal energy functionals so as to close the a priori estimates. In particular, the energy functional for the second component f 2 in the Caflisch's decomposition is complicated, because there is a subtle interplay with f 1 . For this, we make use of the linear combination of estimates for the two functionals, where the smallness of the shear rate α and finiteness of the domain play an important role. Specifically, we obtain estimates (7.2) and (7.3) for the weighted L ∞ norms. To treat L 2 estimates on the right hand side of (7.3), we construct another functional E int (t) in Lemma 7.2, see (7.29), to capture the macroscopic dissipation, and conclude the desired estimates (7.33) and hence (7.36).
Finally, we remark that there have been extensive studies on the stability of shear flow in the multi-dimensional space domain in the context of fluid dynamic equations, cf. [32], in particular, we mention important contributions to the mathematical theories in [6][7][8] by Bedrossian et. al. for either an infinite 2D channel domain T x × R y or an infinite 3D channel domain T x × R y × T z , and an interesting work by Ionescu-Jia [26] for the asymptotic stability of the Couette flow for the 2D Euler equations in the 2D finite channel domain T x × [0, 1] with the zero normal velocities at two boundary planes y = 0, 1, see also the nice survey [5] and references therein. In fact, in comparison with the 1D problem (1.8) under consideration, it would be more interesting to study the existence and asymptotic stability of stationary solutions in the multi-dimensional setting corresponding to those works on fluid dynamic equations. Moreover, it is also challenging to study the fluid dynamic limit for these problems as in [13,18,19] when the vanishing Knudsen number is taken into account. We expect that this paper together with [14] can shed some light on the future investigation on the above problems.
The rest of this paper is organized as follows. In Section 2, we give some basic estimates on the linearized and nonlinear collision operators. In particular, we obtain Lemma 2.4 which is crucially used to obtain the smallness of the nonlocal operator K for large velocity. In Section 3, we revisit Ukai's trace theorem in both the steady and time-dependent cases for the transport operator with the shear force in the 1D setting under consideration. In Section 4 and Section 5, we establish estimates on the first order correction G 1 and the remainder G R , respectively, and hence complete the proof of Theorem 1.1 without showing nonnegativity of the stationary solution. Then we study the time-dependent problem for the local in time existence in Section 6 and the exponential time asymptotic stability of the stationary solution in Section 7 so that the nonnegativity of stationary solution follows. The appendix Section 8 includes some estimates on the boundary product measure when there are multiple bounces induced by the diffuse boundary condition.
Notations. We list some notations and norms used in the paper. Throughout this paper, C denotes some generic positive (generally large) constant and λ denote some generic positive (generally small) constant. D E means that there is a generic constant C > 0 such that D ≤ CE. D ∼ E means D E and E D. 1 A indicates the characteristic function on the set A. We denote · the L 2 ((−1, 1) × R 3 )−norm or the L 2 (−1, 1)−norm or L 2 (R 3 )−norm. Sometimes without any confusion, we use · L ∞ to denote either the L ∞ ([−1, 1] × R 3 )−norm or the L ∞ (R 3 )−norm. Moreover, (·, ·) denotes the L 2 inner product in (−1, 1) × R 3 with the L 2 norm · and · denotes the L 2 inner product in R 3 v . We denote by γ , v y > 0} the incoming set, and by γ 0 = {(±1, v)|v ∈ R 3 , v y = 0} the grazing set. Furthermore |f | 2,± = |f 1 γ± | 2 represent the L 2 norm of f (y, v) at the boundary y = ±1. Finally, we define where n(±1) = (0, ±1, 0). One sees that P γ f defined on {±1} × R 3 , is an L 2 v -projection with respect to the measure |v y | µ(v)dv for any function f defined on γ + .

Basic estimates
In this section we summarize some basic estimates to be used in the following sections. Let us first give some elementary estimates for the linearized collision operator L and nonlinear collision operator Γ, defined by and respectively. It is known that where Q gain denotes the positive part of Q in (1.4). Note that ν 0 is a positive constant in the case of Maxwell molecule collision. The kernel of L, denoted as ker L, is a five-dimensional space spanned by . Define a projection from L 2 to ker L by for g ∈ L 2 , and correspondingly denote the operator P 1 by P 1 g = g − P 0 g, which is orthogonal to P 0 in L 2 . It is also convenient to define with q ≥ 0, then it also holds for any ε ≥ 0 small enough.
For the velocity weighted derivative-in-v x estimates on the nonlinear operator Γ, we have the following lemma.
Lemma 2.2. In the Maxwell molecular case, it holds that

6)
and for any integers m ≥ 0 and q ≥ 0.
Then, by taking directly the L ∞ norm, (2.7) holds because for any integer m ≥ 0. This completes the proof of Lemma 2.2.
The following lemma can be found in [ (2.8) Moreover, for any integer m > 0, there are constants δ 1 > 0 and C > 0 such that Proof. Since (2.8) is quite elementary, we only show (2.9). As in Lemma 2.2, the key point here is to show that the action of the derivatives ∂ m vx on the nonlocal operator L does not involve any other partial derivatives such as ∂ vy or ∂ vz . By (2.1) and (2.8), we have where we have used the change of variableũ = v * − v again. Consequently, as for (2.6), it follows for small enough constants η > 0 and η 1 > 0, where Sobolev's interpolation inequality ∂ m1 vx f 2 ≤ η 1 ∂ m vx f 2 + C η1 f 2 has been used. One the other hand, it can be easily checked that Finally, plugging (2.11) and (2.12) into (2.10) gives (2.9). This completes the proof of Lemma 2.3.
Next, the following lemma which was proved in [14, Proposition 3.1, pp.13] plays a significant role in obtaining the L ∞ estimates of the first component in the Caflisch's decomposition of solutions.
Lemma 2.4. Let K be given by (2.4), then for any nonnegative integer m ≥ 0, there is C > 0 such that for any arbitrarily large q > 0 we have Proof. Since the general case Then, by the similar calculation for estimating I 1 and I 2 in [14, Proposition 3.1, pp.13] yields (2.13). This completes the proof of Lemma 2.4.

A trace theorem
In this section, we present the following Ukai's trace theorem, see also Then, there exists a constant C ε,h > 0 depending on ε and h such that Moreover, it also holds Proof. To prove (3.1), we only consider the case that the boundary phase is outgoing, because the incoming case can be treated similarly. We introduce a parameter t ∈ R and treat (y, v) as Then it follows for (y, v) ∈ γ + \γ ε + . Along this trajectory, one has the identity Next, in light of the Jacobian and by a change of variable Consequently, the desired estimate (3.1) in the case of outgoing boundary follows from (3.8), (3.9) and (3.6).
We now turn to prove (3.2).
Actually, given (t, y, u) for u · n(y f ) > 0. It is easy to see that 0 ≥ s ≥ −t b (y f , u), and it is natural to require that t − s ≤ T. By a change of variable (y, v) → (s, u) and using (3.7), one has On the other hand, if we denotet = t − s, then it follows s ≥ T 1 −t due to t ≥ T 1 . In summary, one has Therefore, we have by change of variable t →t that |f (t + s, Y (t + s;t, y f , u), V (t + s;t, y f , u))|u y |dsdudt. (3.12) Consequently, (3.11) and (3.12) imply (3.10). In addition, it follows that For any (t, y f , u) ∈ [ε 1 , T ] × γ + \γ ε + with ε 1 > 0 to be determined later and for 0 ≥ s ≥ max{−t b (y f , u), ε 1 − t}, we then get from (3.13) and (3.10) that |f (t, y, u))|dtdydu where we have used the fact that hε ≤ t b (y f , u) ≤ h ε due to (y f , u) ∈ γ + \γ ε + . Next, applying Fubini's Theorem and using (3.10) again, one also has Once (3.14) and (3.15) are obtained, it remains now to compute |f (t, y f , u)||u y |dudt.
In fact, if we choose ε 1 to be small enough so that ε 1 ≤ hε, at this stage, the backward trajectory hits the initial plane first. Therefore, for (t, y f , u) ∈ [0, ε 1 ] × γ + \γ ε + , by directly using (3.7) and applying (3.10) once again, it follows The proof of Lemma 3.1 is then completed.

Steady problem: the first order correction
In this and the next sections, we are going to show Theorem 1.1 for the existence of solutions to the steady problem (1.8). Recall (1.9) and (1.10). For the purpose, we will first study in this section the first order correction term G 1 determined by the boundary value problem (1.11) and (1.16). Notice that (1.15) and (1.12) are satisfied. Existence of the remainder G R for the boundary value problem (1.13) and (1.14) will be considered in the next section. Indeed, we have the following proposition. and for any integers m ≥ 0 and q ≥ 0, whereC 1 > 0 is a constant depending only on m and q.
To prove this proposition, let 0 < ǫ < 1 and 0 ≤ σ ≤ 1, then we consider the following general approximation equations and Recall that ν 0 and K are defined by (2.3). The above boundary value problem can be formally reduced to as σ → 1 − and ǫ → 0 + . To prove this rigorously, we deduce the following a priori estimate.
where N 0 is an arbitrary non-negative integer and the constant C 0 > 0 is independent of ǫ and σ.
Proof. The proof of (4.5) is divided into two steps. and We write the solution of (4.6) and (4.7) in the following mild form and we see that

Consequently, we have
where k w is given in Lemma 2.1. Then we iterate (4.10) once more to obtain with By using (4.9) and Lemma 2.1, we see that the last two terms can be bounded as For the other four terms, we only compute I 1,2 because the other three terms can be treated similarly. The estimates are divided into three cases. First of all, we take M > 0 large enough. Case 1. |v| > M . In this case, Lemma 2.1 and (4.9) directly give so that one of the following two estimates holds correspondingly .
This together with Lemma 2.1 and (4.9) gives In this situation, we make use of the boundedness of the operator K on the complement of a singular set. For any large N > 0, we choose a number M (N ) to define The first two difference terms lead to the small contribution of I 1,2 as For the last term, we use the following decomposition which together with Lemma 2.1 as well as (4.9) implies As to I II 1,2 , since y ′′ − y ′ ≤ η 0 , we have that for β ∈ (0, 1), where we have used the fact that Plugging (4.13) into I II 1,2 , we get As a consequence, one has Substituting the above estimates into (4.11), we conclude (4.14) A linear combination of (4.14) from m = 0 to m = N 0 gives the following a priori estimate where C > 0 depends on N 0 and q. This concludes the L ∞ estimate.
L 2 estimates. To close the L ∞ estimate (4.15), we need to derive the L 2 estimate for G 1 . For this, we first consider the zero-th order L 2 estimate G 1 . Notice that G 1 = P 0 G 1 + P 1 G 1 and On the other hand, from (4.
To obtain the L 2 estimate of the macroscopic component, it remains now to deduce the L 2 estimate of b 1,1 . Actually, one can show that where C > 0 is a constant independent of ǫ and σ. For this, we define Taking the inner product of (4.3) and Ψ b1,1 over (−1, 1) × R 3 , we get We now compute the terms in (4.19) one by one. By Cauchy-Schwarz's inequality and (4.18), one has For the boundary term, one has from (4.4) and (4.18) that Combining the above estimates for the terms in (4.19), we have (4.17). We now deduce the L 2 estimate on the microscopic component P 1 G 1 . Direct energy estimate on (4.3) gives Thus, (4.17) and (4.20) as well as (4.16) yield Furthermore, for the higher order L 2 estimates on where Lemma 2.3 has been used for σ(∂ m vx LG 1 , ∂ m vx G 1 ). Finally, the a priori estimate (4.5) follows from (4.15), (4.21) and (4.22). This completes the proof of Lemma 4.1.
With the a priori estimate (4.5), we now prove the following existence result for general linear equations (4.3) and (4.4). Before doing this, we first define the following function space And for convenience, we also define a linear operator L σ by then there exists a unique solution G 1 = G 1 (y, v) to (4.3) and (4.4) with σ = 1 satisfying where C 0 > 0 is a constant depending only on N 0 and q.
Proof. The proof is based on a bootstrap argument in the following three steps.
Step 1. Existence for σ = 0. If σ = 0, then (4.3) and (4.4) are reduced to which has a unique explicit solution (4.23), and a direct calculation implies Step 2. Existence for any σ ∈ [0, σ * ] with some σ * > 0. Suppose σ ∈ (0, 1], we consider a more general equation To solve this boundary value problem, we design the following approximation equations is well defined by step 1 and satisfies the estimate for any n ≥ 0. Furthermore, one can also show that for σ ∈ [0, σ * ], there exists a unique solution G 1 ∈ X N0 to the problem (4.26) and (4.27). Actually, the a priori estimate (4.5) implies that we still have the bound In other words, it follows Step 3. Existence for σ ∈ [0, 2σ * ]. By using (4.31) and performing similar calculations as for obtaining (4.29) and (4.30), one can see that there exists a unique solution G 1 ∈ X N0 to the lifted equation . Therefore, the solution mapping L −1 2σ * is also well-defined on X N0 and the estimate (4.24) holds for σ = 2σ * .
Finally, repeating the above procedure step by step, one can reach σ = 1 so that L −1 1 exists and (4.24) also follows simultaneously. This completes the proof of Lemma 4.2.
Proof of Proposition 4.1: By setting F = −v x v y √ µ in Lemma 4.2, we see that for any ǫ > 0 there exists a unique solution G ǫ 1 ∈ X N0 to the boundary value problem 1 satisfies (4.1) and the estimate G ǫ 1 XN 0 ≤C 1 , whereC 1 > 0 is independent of ǫ. Furthermore, we define a positive sequence {ǫ n } ∞ n=1 such that |ǫ n+1 − ǫ n | ≤ 2 −n , then ǫ n → 0 + as n → +∞. We consider the following approximation equations is a Cauchy sequence in X N0 . Thus, letting n → ∞, the limit function denoted by G 1 is the unique solution of (1.11) and (1.16). Moreover, G 1 satisfies (4.1) and the bound (4.2). The proof of Proposition 4.1 is then completed.

Steady problem: remainder
Based on Proposition 4.1, one can further show the following existence result for the remainder G R . Recall the steady problem (1.8) as well as (1.9) and (1.10).
Proposition 5.1. The boundary value problem (1.13) and (1.14) admits a unique solution And there is an integer q 0 > 0 such that for any integer q ≥ q 0 , there is α 0 = α 0 (q) > 0 depending on q such that for any α ∈ (0, α 0 ) and any integer m ≥ 0, it holds that w q ∂ m vx G R L ∞ ≤C R , whereC R > 0 is a constant depending only on m and q but independent of α.

5.1.
Caflisch's decomposition. To prove Proposition 5.1, we follow the strategy of the proof in [14] for treating the shear force term in the framework of perturbation. In fact, notice that there is a growth term α 2 v x v y G R in the equation (1.13). To treat this growth in velocity, the key point is to use the Caflisch's decomposition [10] and an algebraic weighted estimate introduced originally by Arkeryd-Esposito-Pulvirenti [4]. For the purpose, we first decompose the remainder where G R,1 and G R,2 satisfy the following two boundary value problems, respectively, and Here χ M (v) is a non-negative smooth cutoff function such that and K is defined by (2.4). The existence of (5.2), (5.3), (5.4) and (5.5) can be constructed via the approximation sequence by iteratively solving the following system and for a small parameter ǫ > 0, where we have set [G 0 R,1 , G 0 R,2 ] = [0, 0] for n = 0. The proof of Proposition 5.1 follows by three steps. First, similarly for treating the existence of G 1 , we introduce a modified coupled boundary value problem with two parameters ǫ > 0 and 0 ≤ σ ≤ 1. This boundary value problem is directly solvable via the characteristic method in case of σ = 0 corresponding to the homogeneous inflow data, and we then lift the value of σ from σ = 0 for the zero inflow data to σ = 1 for the full diffuse reflection boundary condition by a bootstrap argument. Second, we establish the limit n → +∞ for any fixed parameter ǫ > 0. Third, we pass the limit ǫ → 0 + to obtain the desired solution which satisfies (5.2), (5.3), (5.4) and (5.5). As a result, with the help of (5.1), we get the solution to the original boundary value problem (1.13) and (1.14).

5.2.
A priori estimates with parameters ǫ and σ. Let us first show that [G n+1 R,1 , G n+1 R,2 ] is well-defined once [G n R,1 , G n R,2 ] is given. To do this, we apply the method of contraction mapping. We define the linear vector operator parameterized by σ ∈ [0, 1] as follows: We then consider the solvability of the following coupled linear system where F 1 , F 1 and F 2,b are given, and F 1 , 1 + F 2 , √ µ = 0. In the rest of the proof, for brevity, we denoteF In what follows, we look for solutions to the system (5.10) in the Banach space supplemented with the norm Let us now deduce the a priori estimate for the parameterized linear system (5.10).
Proof. The proof is divided into two steps.
where n(y j ) = (0, 1, 0) if y j = 1 and n(y j ) = (0, −1, 0) if y j = −1. Let the iterated integral for k ≥ 2 be defined as where dσ l = √ 2πµ(v l )|v ly |dv l is a probability measure. Without lost of generality, we assume lim provided that ǫ > 0 and qα > 0 are suitably small. By Lemma 2.4, it is straightforward to see By taking q sufficiently large, (5.21) further gives Similarly, one can also write the solution of (5.15) and (5.16) in the mild form of , and for k ≥ 2, Vj (w q F 2,b )(t l , y l , V l−1 cl (t l )) dΣ l (t l ), Moreover, in the above expressions we have used the following notations A ǫ (τ,V l cl (τ ))dτ dσ j , The L ∞ estimates for H 2,m is more complicated because K has no smallness property. To overcome this, we have to iterate (5.23) twice. Let us first compute I n (1 ≤ n ≤ 10) term by term. Recalling the definition (2.5) for k w , one directly has by (5.20) By Lemma 2.2, it follows and similarly, It is straightforward to see In the sequel, for simplicity, we denote C m,q for the constant m!4 q q!. By Lemma 8.1, it follows where we have taken T 0 = t − t k with k = C 1 T 5 4 0 , and both C 1 > 0 and C 2 > 0 are given in Lemma 8.1.
Putting all the estimates for I n (1 ≤ n ≤ 10) above together and adjusting the constants, we have where Then let us define a new backward time cycle as and the starting point , for some s ∈ R and l ∈ Z + . Furthermore, for ℓ ∈ Z + , we also denotē . Iterating (5.26) again, one has where according to Lemma 8.1 and Remark 8.1, we choose ı ∈ Z + such that ı ∼ (T 0 ) 5 4 with T 0 = s − t ′ ı being suitably large. We claim that where η > 0 is suitably small. To prove (5.28), we only estimate the fourth term on the right hand side of (5.27), because the other terms can be estimated similarly. For any sufficiently small , then rewrite the fourth term on the right hand side of (5.27) as By Lemma 8.1, it is easy to see For J 1 , the computation is divided into the following three cases.
In this case, by Lemma 2.1, it follows that where C(q) > 0 and depends on q!. Therefore, one has by using Lemma 8.1 again Note that here and in the sequel, l and ℓ run over [1, k − 1] and [1, ı − 1], respectively.
Then either of the following two estimates holds correspondingly .
This together with Lemma 2.1 imply The key point here is to convert the L 1 integral with respect to the double v variables into the L 2 norm with respect to the variables y and v. To do so, for any large N > 0, we choose a number M (N ) to define k w,M (u, v ′ ) as (4.12), then decompose . From Lemma 8.1, the first two difference terms lead to a small contribution in J 1 as For the remaining main contribution of the bounded product k w, and apply a change of variable v ′ ly →ỹ. Then one has We now estimate this part as follows Putting all the estimates for J 1 and J 2 together, we now obtain As mentioned before, by performing the similar calculations for the other terms on the right hand side of (5.27), one has Since ı ∼ (T 0 ) 5 4 , by taking M and N large enough and η 0 = (T 0 ) − 5 2 small enough, (5.29) further yields (5.28). Finally, taking a linear combination of (5.28) with m = 0, 1, · · · , N 0 , we conclude Remark 5.1. We point out that the estimates (5.22) and (5.30) obtained above are independent of ǫ. Moreover, bothT 0 and T 0 are independent of t, because starting from any t ∈ (−∞, +∞) we can trace back k times to some t k which can be negative.
Step 2. L 2 estimate. To close the final estimate, we turn to obtain the L 2 estimate of ∂ m vx G 2 with 0 ≤ m ≤ N 0 . The goal is to prove that for given ǫ > 0 there exists C(ǫ) > 0 depending on ǫ such that For this, we begin with the following equations for G 2 Taking the inner product of (5.32) 1 and G 2 over (y, v) ∈ (−1, 1) × R 3 , we have, for η > 0 where the following estimate on the boundary term has been used Here we have used the notation and the estimate vy ≷0 G 1 (±1)|v y |dv ≤ C w q G 1 L ∞ , for q > 5/2.

Next, since
by dividing the domain for integration as one sees that the grazing part of |P γ G 2 (1)| 2 2,+ is bounded by the Hölder inequality as For non-grazing region, we have by using the trace Lemma 3.1 that Putting (5.34) and (5.35) together, one has Consequently,(5.33) and (5.36) give Remark 5.2. Note that the constant C(ǫ) in (5.37) is independent of the parameter σ.
It remains to deduce the L 2 estimate of the higher order velocity derivatives. For this, applying ∂ m vx (m ≥ 1) to (5.32) Taking the inner product of (5.38) and ∂ m vx G 2 , we deduce where we have used the following estimates for the incoming boundary term by (5.38 Proof. The proof is based on the a priori estimate (5.12) established in Lemma 5.1 and a bootstrap argument. As for Lemma 4.2, the proof is also divided into the following three steps.
Step 1. Existence for σ = 0. If σ = 0, then (5.10) is reduced to respectively. Then, in this simple case, the existence of L ∞ -solutions can be directly proved by the characteristic method to have Step 2. Existence for σ ∈ [0, σ * ] for some σ * > 0. Let σ ∈ (0, 1], we now consider and For the above system, we design the following approximation scheme The goal in the following proof has twofold: whereC 1 > 0 is independent of σ and n. Choosing 0 < σ * < 1 suitably small such that which is equivalent to Step 3. Existence for σ ∈ [0, 2σ * ] for some σ * > 0. By using (5.54) and performing the similar calculations as for obtaining (5.52) and (5.53), for σ ∈ [0, σ * ], one can see that there exists a unique solution [G 1 , G 2 ] ∈ X α,N0 to the lifted system and In other words, we have shown the existence of L −1 2σ * on X α,N0 and (5.12) holds true for σ = 2σ * . Therefore, by repeating this procedure in finite time, , one can see that L −1 1 exists when σ = 1 and (5.40) follows correspondingly. This completes the proof of the lemma. Let us first go back to the approximation system (5.6), (5.7), (5.8) and (5.9). By applying Lemma 5.2, for fixed ǫ > 0, we see that [G n+1 R,1 , G n+1 R,2 ] is well defined when [G n R,1 , G n R,2 ] is given and the solution belongs to X α,N0 defined in (5.11) for N 0 ≥ 0.
We now show that {G n R,1 , G n R,2 } ∞ n=0 is a Cauchy sequence in X α,N0 , which hence implies that its limit denoted by [G ǫ R,1 , G ǫ R,2 ] is the unique solution of the following system and Furthermore, we will show that the convergence of the sequence {G n R,1 , G n R,2 } ∞ n=0 is independent of ǫ. For this, we first prove that where C 0 > 0 is independent of ǫ and n for all n ≥ 0. We apply induction in n. Notice [G 0 R,1 , G 0 R,2 ] = [0, 0]. If n = 1, then the system (5.6), (5.7), (5.8) and (5.9) reads and where the constant C > 0 is independent of ǫ, see also Remark 5.1.
Since the mass of G 1 R,2 is not conserved, in order to estimate the macroscopic component of G 1 R,2 we instead turn to obtain the L 2 estimate of G 1 By (5.60), (5.61), (5.62) and (5.63), it is easy to see that G 1 R satisfies and Next, for n ≥ 1, denote and define the projectionP 0 , from L 2 to ker L, as In addition, we will also use the following notations . Note that a n = a n 1 + a n 2 , b n = b n 1 + b n 2 , c n = c n 1 + c n 2 , 1 −1 a n (y)dy = 0.
for q > 5/2, to obtain the estimate of [a 1 2 , b 1 2 , c 1 2 ] , it suffices to derive the L 2 estimates of [ a 1 , b 1 , c 1 ]. In what follows, we will show that the L 2 norm of the macroscopic part of G 1 R can be indeed dominated by its microscopic component and other known terms. We estimate [a 1 , b 1 , c 1 ] by the dual argument. First of all, we let Ψ(y, v) ∈ C ∞ ([−1, 1] × R 3 ), and take the inner product of (5.66) and Ψ over (−1, 1) × R 3 , to obtain Estimate on a 1 . Let Plugging Ψ = Ψ a 1 into (5.71), we now compute the equation term by term. First of all, by Cauchy-Schwarz inequality with η > 0 and using (5.73), one has Then by Lemmas 4.1 and 2.2, it follows has by a similar argument as above that The last boundary term v y G 1 R (1), Ψ a 1 (1) − v y G 1 R (−1), Ψ a 1 (−1) vanishes because of (5.72). Putting all the estimates above together, we have We now compute each term in (5.71) with Ψ = Ψ b 1 i . By Cauchy-Schwarz inequality and (5.77), one has Similar to (5.74) and (5.75), it follows For the boundary term, noting that , Ψ bi (±1) = 0 has been used. We now conclude from the above estimates for b 1 Estimate on c 1 . Let By Cauchy-Schwarz inequality and (5.80), it follows Also similar to derive 5.74 and (5.75), one has (µ − 1 2 S 0 , Ψ c 1 ) ≤ η c 1 2 + C, and For the boundary term, by applying (5.78) and using Combining the above estimates on c 1 give Finally, a linear combination of (5.76), (5.79) and (5.81) gives This together with (5.69) and (5.70) implies that [a 1 In order to obtain the estimates for P 1 G 1 R,2 , we have to further consider the BVP for G 1 R,2 as follows: Applying the estimates (5.33), (5.36) and (5.39) with σ = 1, F 2 = 0 and F 2,b = 0, one has and where all the constants on the right hand side are independent of ǫ. Then (5.82), (5.84), (5.85) and (5.86) give Consequently, a linear combination of (5.64), (5.65) and (5.87) gives for some suitably large C 0 > 0. Therefore (5.59) holds for n = 1. We now assume that (5.59) is valid for n = k ≥ 1 and then prove that it holds for n = k + 1. In fact, applying the estimates (5.22) and (5.30) to the system (5.6)-(5.7) and (5.8)-(5.9) with n = k, one has . By Lemma 4.1, Lemma 2.2 and the induction assumption, we have For the L 2 estimate, by performing a parallel calculation as for (5.83), one has and Ψ c k+1 in the same way as for Ψ a 1 , Ψ b 1 i and Ψ c 1 , respectively. Hence, (5.91) also gives and provided that α is chosen to be sufficiently small. Thus (5.59) holds for n = k + 1. Therefore, (5.59) holds for all n ≥ 0. We now turn to prove that [G n R,1 , G n R,1 ] ∞ n=0 is a Cauchy sequence in X α,N0 . For this, denote We claim that On the other hand, we have from Lemma 2.2 that where C > 0 is independent of ǫ. Furthermore, by taking the limit ǫ → 0, we can repeat the same procedure like letting n → ∞ so that the limit function [G R,1 , G R,2 ] ∈ X α,N0 is the unique solution of (5.2)-(5.3) and (5.4)-(5.5) with the same bound as (5.98). Then the proof of Proposition 5.1 is completed.
Finally, Theorem 1.1 is an immediate consequence of Proposition 4.1 and Proposition 5.1, except for the non-negativity of the solution F st (y, v) that will be proved from the dynamical stability of F st (y, v) in Theorem 1.2.

Unsteady problem: local existence
We now turn to the time-dependent situation. To solve the initial boundary value problem (1.21), we set the perturbation as The goal of this section is to construct the local-in-time solution to the initial boundary value problem (6.2). The proof of the global existence of solutions as well as the large time behavior will be left to the next section. To resolve the difficulty caused by the growth term α 2 v x v y f , it is still necessary to introduce the decomposition where f 1 and f 2 satisfy the following initial boundary value problem and respectively. Note that initial data for f 2 is set to be zero. We will look for solutions to (6.4)-(6.5) and (6.6)-(6.7) in the following function space associated with the norm Proof. We first consider the following system for approximation solutions and where for any m ≥ 0, provided that This also implies that [f n+1 ] is well-defined by (6.8)-(6.9) and (6.10)-(6.11) if [f n 1 , f n 2 ] is bounded as in (6.12). Denote , √ µG n = G n 1 + √ µG n 2 , then G n 1 and G n 2 satisfy and Here, andw 2 is given by (5.25).
For the boundary terms I and I (2) b , we use the equations (6.17) and (6.18) recursively to obtain l (s)ds, l (s)ds, l (s)ds, and I 3 =W where k ≥ 2. Here, similar to (5.24),Σ (i) l (s) (i = 1, 2) is given as To obtain (6.12), one can first prove that for fixed finite k > 0 and any t ≥ 0, by choosing C 0 > 0 suitably large. Note that (6.19) can be easily obtained by using (6.17) and (6.18) recursively because k is finite.
In the following, we prove (6.12) for m = n + 1 under the assumption that it holds for m ≤ n. By letting t ≤ T * with T * > 0 being suitably small and applying Lemma 8.1, we have and sup 0≤t≤T * where Lemma 2.4 has been used to have the factor 1 q in (6.20), and the coefficient T * on the right hand side of (6.21) comes from the last two terms in (6.18) as well as I 3 . Choosing T * andε suitably small so that C(T * +ε) ≤ 1 8 , and using the induction argument, we have from (6.21) and (6.20) that where [ n k ] stands for the largest integer no more than n k . Therefore, , and √ µG n =G n 1 + √ µG n 2 . Then similar to (6.20) and (6.21), one has By taking q > 0 sufficiently large and α > 0, ε 0 > as well T * > 0 suitably small, we have from (6.24) and (6.22) that On the other hand, sup 0≤l≤k sup 0≤t≤T * [G l 1 ,G l 2 ] L ∞ is bounded due to (6.12). Hence it follows that {[f n 1 , f n 2 ]} ∞ n=1 is a Cauchy sequence in the space Y α,T * , and there is a unique [f 1 , f 2 ] ∈ Y α,T * such that [f n 1 , f n 2 ] converges strongly to [f 1 , f 2 ] as n → +∞ and [f 1 , f 2 ] is the desired local in time solution to the coupled system (6.4)-(6.5) and (6.6)-(6.7). This completes the proof of Theorem 6.1.

Unsteady problem: asymptotic stability and positivity
This section is about the global existence and large time behavior of solution to the initial boundary value problem (6.2). Recall the decomposition (6.3) with f 1 and f 2 satisfying the coupled system (6.4)-(6.5) and (6.6)-(6.7). Firstly, we focus on the uniform L ∞ ∩ L 2 estimates under the a priori assumption sup s≥0 {e λ0s w q f 1 (s, y, v) L ∞ + e λ0s w q f 2 (s, y, v) L ∞ } ≤ε, (7.1) for a constantε > 0 suitably small, where λ 0 > 0 independent of α is to be determined later. And then we will give the proof of Theorem 1.2.
7.1. L ∞ estimates. As in the proof of Theorem 6.1, the L ∞ estimates of f follows from the uniform L ∞ estimates on f 1 and f 2 .
Proof. For brevity, set [g 1 , g 2 ](t, y, v) = e λ0t w q (v)[f 1 , f 2 ](t, y, v) (7.4) with λ 0 > 0 to be chosen. Then, the IBVP for [g 1 , g 2 ] is given as follows: Along the characteristic line (3.3) the solution to the above problem can be written in the mild form: and where Consequently, for any t ≥ 0, by applying Lemmas 2.2, 2.4 and 8.1 as well as the a priori assumption (7.1), we get from (7.5) that that gives (7.2). For g 2 , similar to (5.26), one has l (s)ds + P(t), where We now have by iterating (7.7) that Next, combining (7.2) at t = T 0 and (7.10), one has Then it follows that for any t ∈ [0, T 0 ], In particular, we have Moreover, (7.11) can be extended to for any t ∈ [s, s + T 0 ] with s ≥ 0.
Next, for any integer m ≥ 1, we can repeat the estimate (7.12) in finite times so that the functions [f 1 , f 2 ](lT 0 + s) for l = m − 1, m − 2, ..., 0 satisfy 0≤s≤mT0 e λ0s f 2 (s) . (7.14) Furthermore, for any t ≥ T 0 , we can find an integer m ≥ 0 such that t = mT 0 + s with 0 ≤ s ≤ T 0 . Then we have, on one hand, by (7.14), that On the other hand, (7.13) implies that which is equivalent to Consequently, applying (7.15) to (7.16) gives the second estimate (7.3). This together with (7.2) concludes the L ∞ estimate on f 1 and f 2 , and then it completes the proof of Lemma 7.1.
By putting the above estimate back to (7.3) and using the smallness of α andε, one has sup 0≤s≤t e λ0s w q f 2 (s) L ∞ ≤ C w q f 0 L ∞ + C sup 0≤s≤t e λ0s w q f 1 (s) L ∞ . (7.36) Moreover, by plugging (7.36) to (7.2) and using the smallness of α andε as well as (7.36), it holds that sup 0≤s≤t e λ0s w q [f 1 , f 2 ](s) L ∞ ≤ C w q f 0 L ∞ , which gives (1.23). Since w q f 0 L ∞ is sufficiently small, the a priori assumption (7.1) is closed. Finally, the non-negativity of the global solution constructed above can be proved similar to [14] so that the proof of Theorem 1.2 is completed.

Appendix
Recall the backward time cycle starting at (t 0 , y 0 , v 0 ) = (t, y, v) in (5.18), the boundary probability measure dσ l on V l in (5.19) and the product measure dΣ l (s) over k−1 j=1 V j in (5.24). The following lemma gives an estimate on the measure of the phase space Π k−1 j=1 V j when there are k times bounce. In particular, let T 0 > 0 large enough, there exist constants C 1 and C 2 > 0 independent of T 0 such that for k = C 1 T Here we have used by the Peetre's inequality and the fact that |V j cl (t j+1 ) − v j | = α|t b (v j )v jy | ≤ 2α. Then the proof of lemma is completed.