{'fraction_non_alphanumeric': 0.05803511292712496, 'fraction_numerical': 0.023214045170849983, 'mean_word_length': 3.8180638182801516, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A convergent method for linear half-space kinetic equations

We give a unified proof for the well-posedness of a class of linear half-space equations with general incoming data and construct a Galerkin method to numerically resolve this type of equations in a systematic way. Our main strategy in both analysis and numerics includes three steps: adding damping terms to the original half-space equation, using an inf-sup argument and even-odd decomposition to establish the well-posedness of the damped equation, and then recovering solutions to the original half-space equation. The proposed numerical methods for the damped equation is shown to be quasi-optimal and the numerical error of approximations to the original equation is controlled by that of the damped equation. This efficient solution to the half-space problem is useful for kinetic-fluid coupling simulations.


Introduction
In this paper we propose a Galerkin method for computing a class of half-space kinetic equation with given incoming data: Here u ∈ R is a given constant and V ⊆ R ν with ν ≥ 1. We are mainly concerned with several cases for V: V = [−1, 1], V = S ν−1 or V = R ν . Typical examples are V = [−1, 1] or V = S 2 (with u = 0) for neutron transport equations, V = R 3 (with arbitrary u ∈ R) for linearized Boltzmann equations, and V = R 1 , R 3 for linearized BGK equations. If V = S 2 , then in spherical coordinates v = (cos φ, cos θ sin φ, sin θ sin φ) ∈ S 2 and v 1 = cos φ. If V = R 3 , then v = (v 1 , v 2 , v 3 ) ∈ R 3 and v 1 is the first component of v. Our framework also covers systems with multiple species, in which case, the ambient space for the density function f is (L 2 ( dv)) d with d ≥ 1, where d is the number of species and dv is the usual Lebesgue measure for R d or the rotation-invariant measure on S d−1 . The linear operator L only acts on the velocity variable v and is nonlocal in v in general. The specific structure and main assumptions regarding L will be given in Section 2. While the half-space kinetic problem is an important problem on its own, the study of such equations in the form of (1.1) is also our first step in developing systematic algorithms that resolve couplings between kinetic and moment-closure equations. This type of couplings are common in domain-decomposition methods for solving kinetic equations with multi-scales. Multi-scale phenomenon naturally arises in gas dynamics where dense and dilute fluids coexist in a region. It is also common in transport phenomenon where the background medium has different scattering and absorption properties. Taking gas dynamics as an example, in a typical domain-decomposition method, one uses the original kinetic equation for the dilute part and a fluid equation for the dense side. The boundary conditions for the two sides are encoded in the density function near their interface. Due to the sharp variation in scales, the density function in general will have a finite jump across the interface, thus boundary layers will form there. Half-space equations are exactly leading-order boundarylayer equations whose solutions bridge the gap between the fluid and kinetic boundary conditions. Therefore understanding the well-posedness of (1.1) and constructing accurate and efficient numerical schemes to resolve this equation will provide explicit characterization of the couplings. We note that due to the lack of efficient and accurate solvers of the half-space equation, various methods have been developed to bypass solving the half-space equation (see for examples [DJM05,LM12]). Our goal instead is to resolve the kineticfluid coupling through efficient and accurate solutions of the half-space equations, which will be the focus of our forthcoming work.
In the literature the well-posedness of equation (1.1) has been referred to as the Milne problem and has long been investigated [BSS84,BY12,CGS88,Gol08,UYY03] for various models. For example, for the conservative neutron transport equation where the null space of L consists of only constant functions, it has been shown [BSS84,BLP79] that there exists a unique solution f ∈ L ∞ ( dx dv) and the end-state of f at x = ∞ is a constant. In general for models where dim(Null L) > 1, additional boundary conditions are needed for (1.1) to have a unique solution [CGS88]. A typical example is the linearized Boltzmann equation of one species over v ∈ R 3 . In this case we have dim(Null L) = 5. The well-posedness of this equation with hard-sphere collisions is fully resolved in the fundamental work by Coron, Golse, and Sulem [CGS88]. In this work, it is shown that depending on the choices of u, one needs to prescribe various numbers of additional boundary conditions such that (1.1) is well-posed. These numbers of boundary conditions correspond to the counting of the incoming Euler characteristics at x = ∞. The proof in [CGS88] relies mainly on the energy method. Subsequently, a different proof using a variational formulation of (1.1) for the linearized Boltzmann equation is given in [UYY03]. The key idea in [UYY03] is to revise (1.1) by adding certain damping terms. The revised collision operator thus obtained is coercive and it enforces the end-state of f at x = ∞ to be zero. By the conservation properties of L, the authors then show that (1.1) is well-posed for a large class of incoming data. One restriction in [UYY03] is that u cannot be chosen in the way such that the Mach number of the system is −1, 1, or 0. This restriction was later removed in [Gol08].
The variational formulation is also a common tool in proving the well-posedness of the neutron transport equations over general bounded domains Ω in R ν x . There is a vast literature in this direction and we will only review some of the main framework and results in [ES12] which are most relevant to us. Specifically, suppose Ω ⊆ R ν x is a bounded domain with a C 1 boundary such that the normal direction on ∂Ω is well-defined. Consider the neutron transport equation with the incoming boundary: where n is the normal direction on ∂Ω and v ∈ S ν−1 is a unit vector. The scattering operator L is given by where c a and c s are the absorption and scattering coefficients and σ is the scattering kernel. In [ES12] it is assumed that 0 < α 0 < c a < α 0 and 0 ≤ c s ≤ c 0 , thus L is subcritical and Null L = ∅. The main novelty of [ES12] is that one decomposes the solution f into its even and odd parts in v and imposes different regularities for these two parts. Using this mixed regularity, the authors of [ES12] write (1.2) into a variational form and verify that the bilinear operator involved satisfies an inf-sup condition over a properly chosen function space. Moreover, they show that for appropriately constructed Galerkin approximations, the bilinear operator satisfies the inf-sup condition over finite-dimensional approximation spaces as well. This then shows the Galerkin approximation is quasi-optimal. We also note that the even-odd parity was widely used for transport equations, see for example [JPT01]. There are two main goals in our paper: First, we will generalize the analysis in [ES12, Gol08, UYY03] to obtain a unified proof for the well-posedness of half-space equations in the form of (1.1). Second, we will develop a systematic Galerkin method to numerically resolve (1.1) and obtain accuracy estimates for our scheme.
We now briefly explain our main results and compare them with previous ones in the literature. In terms of analysis, as in [CGS88], we show that with appropriate additional boundary conditions at x = ∞, equation (1.1) has a unique solution. The basic framework we use is the even-odd variational formulation in [ES12]. Compared with [ES12], here we allow the linear operator L to have a nontrivial null space and the background velocity u to be any arbitrary constant for general models. The number of additional boundary conditions will change with u. Due to the loss of coercivity of L, if one directly applies the variational method in [ES12] then the bilinear operator B ceases to satisfy the inf-sup condition. To overcome this degeneracy, we utilize the ideas in [Gol08,UYY03] by adding damping terms to (1.1) and reconstructing solutions to (1.1) from the damped equation. In the case of linearized Boltzmann equation with a single species, we thus recover the results (in the L 2 spaces) in [Gol08,UYY03]. The main differences between our work and [Gol08,UYY03] are: First, we use a different variational formulation which is convenient for performing numerical analysis. Second, the reconstruction in [Gol08,UYY03] is restricted to a set of incoming data with a finite codimension such that the damping terms are identically zero. Here we use slightly different damping terms and we recover solutions to (1.1) from the damped equation for any incoming data. On the other hand, our main concern is the convergence and accuracy of the numerical scheme and the basic L 2 -spaces are sufficient for this purpose. Therefore, except for the hard sphere case, we do not try to achieve decay estimates of the half-space solution to its end-state at x = ∞, while in the literature there are a lot of works that show subexponential or superpolynomial decay of the solution to its end-state for hard or soft potentials for the linearized Boltzmann equation(see for example [CLY04,WYY06,WYY07]).
Our analysis also applies to linearized Boltzmann equations with multiple species and linear neutron transport equations with critical or subcritical scatterings, thus providing an alternative proof to the wellposedness result (in the L 2 -space) in [BY12].
In parallel with the analysis, numerically we first solve the damped half-space equation and then recover the solution to the original equation. We will use a spectral method and achieve quasi-optimal accuracy (for the damped equation) as in [ES12]. The spectral method dates back to Degond and Mas-Gallic [DMG87] for solving radiative transfer equations, and was later extended by Coron [Cor90] to solving the linearized BGK equation as well. Compared with these works, our approach differs in three ways: First, as a result of using the even-odd formulation, we can derive explicit boundary conditions for the approximate equations. In particular, the number of these boundary conditions is shown to be consistent with the number of the unknowns. Hence our discrete systems are always well-posed. This was not the case in [Cor90] where a least square method was used to solve a potentially overdetermined problem. Second, the method in [Cor90] used Hermite polynomials defined on the whole velocity space as their basis functions. This leads to severe Gibbs phenomenon, since in general the solution to the half-space equation has a finite jump at x = 0 and v = −u. Here we choose to use basis functions with jumps at v = −u which naturally fit into the even-odd formulation. Third, we will treat the cases with arbitrary bulk velocities u in a uniform way while in [Cor90] different schemes are used for the cases u = 0 and u = 0.
There are also non-spectral methods developed for solving the half-space equations. For example, the work by Golse and Klar [GK95] uses Chapman-Enskog approximation with diffusive closures. The accuracy of these approximations would be hard to analyze: the iterative approach couples the error from the systematic expansion truncation with the numerical error. Moreover, this work ( [GK95]) also treats the cases u = 0 and u = 0 separately. The very recent work by Besse et al. [BBG + 11] treats the half-space problem as a boundary layer matching kinetics with the limiting fluid equation, where a Marshak type approximation [Mar47] is applied for boundary fluxes. Similar idea was also used in [Del03]. As shown already in [Cor90], in general the Marshak approximation does not yield accurate approximations to the half-space problem.
The layout of this paper as follows: in Section 2, we gather the basic information related to the linear operator L and set up the variational formulation. In Section 3, we summarize the main results in this paper which include both the analytical results and the numerical scheme. We give proofs in Section 4. Section 5 presents more details of the numerical scheme and numerical results for linearized BGK and linear transport equations.

Basic Setting
In this section, we collect the basic information about the collision operator L and present the variational formulation of a damped version of (1.1).
2.1. Linear collision operator. In order to state the main assumptions imposed on L, we first introduce some notations. Denote Null L as the null space of L. Let P : (L 2 ( dv)) d → Null L be the projection onto Null L. Define the weight function for some 0 ≤ ω 0 ≤ 1. The main assumptions on L are as follows: (A1) L : (L 2 ( dv)) d → (L 2 ( dv)) d is self-adjoint, nonnegative with its domain given by where a(v) is defined in (2.1). (A2) L : (L 2 (a dv)) d → (L 2 ( 1 a dv)) d is bounded, that is, there exists a constant C 0 > 0 such that Null L is finite dimensional and Null L ⊆ (L p ( dv)) d for all p ∈ [1, ∞). (A4) L has a spectral gap: there exists σ 0 > 0 such that where P ⊥ = I − P is the projection onto the null orthogonal space (Null L) ⊥ .
One operator that is of particular importance is P 1 : Null L → Null L which is defined by Note that P 1 is a symmetric operator on the finite dimension space Null L. Therefore, its eigenfunctions form a complete set of basis of Null L. Denote H + , H − , H 0 as the eigenspaces of P 1 corresponding to positive, negative, and zero eigenvalues respectively and denote their dimensions as Let X +,i , X −,j , X 0,k be the associated unit eigenfunctions with 1 ≤ i ≤ ν + , 1 ≤ j ≤ ν − , and 1 ≤ k ≤ ν 0 for ν ± , ν 0 = 0. Note that if any of ν ± , ν 0 is equal to zero, then we simply do not have any eigenfunction associated with the corresponding eigenspace. By their definitions, these eigenfunctions satisfy where α ∈ {+, −, 0}, γ ∈ {i, j, k}, 1 ≤ i ≤ ν + , 1 ≤ j ≤ ν − , and 1 ≤ k ≤ ν 0 . Using the notation of H ± , H 0 , we can now state the full equation that we want to study in this paper. Suppose L is a linear operator in v that satisfies (A1)-(A4). Our goal is to prove the well-posedness of, construct efficient numerical schemes, and obtain accuracy estimates for the following equation: The particular formulation about the end-state f ∞ was given in [CGS88] where the authors proved the well-posedness of the half-space linearized Boltzmann equation: CGS88]). Let L be the linearized Boltzmann operator with a hard-sphere collision kernel. Then there exists a constant β > 0 and a unique f ∞ ∈ H + ⊕ H 0 such that equation (1.1) has a unique solution f which satisfies where a(v) = 1 + |v|.
Remark 2.1. The main result in [CGS88] is actually stronger than Theorem 2.1 where f − f ∞ is shown to be in L ∞ (e 2βx dx; L 2 ( dv)). Here we content ourselves with the L 2 -weighted space (in x) since L 2 suffices our needs in proving the quasi-optimal convergence of our numerical scheme.
Many well-known linear or linearized kinetic models satisfy the above assumptions for the collision operators. These include the classical linearized Boltzmann equations for either single-species system or multispecies with hard-sphere collisions and the linear neutron transport equations. The particular equations that we use as numerical examples are the isotropic neutron transport equation (NTE) with slab geometry and the linearized BGK equation. Similar analysis can be carried out to models satisfying (A1)-(A4) without extra difficulties. The main structure of these two equations are as follows. The linear operator of the isotropic NTE is the simplest scattering operator which has the form In this case, a(v) = 1 + |v| = O(1) and (L 2 (a dv)) d coincides with (L 2 ( dv)) d . The linearized BGK operator is the linearization of the nonlinear BGK operator, which is introduced as a simplified model that captures some fundamental behavior of the nonlinear Boltzmann equation. The collision operator of the nonlinear BGK is defined as where M[F ] is the local Maxwellian associated with F defined by For a given bulk velocity u ∈ R, define the global Maxwellian with the steady state (ρ, u, θ) = (1, u, 1/2) as Linearizing the operator Q around M by setting we obtain the linearized BGK operator: where m u (v) is f projected onto the kernel space of L u . In the case of the 1D linearized BGK, one has: Therefore, m u (v) is a quadratic function associated with a Maxwellian to 1/2 power: where ( ρ, u, θ) are defined in the way such that first three moments of m(v) agree with those of f : The half-space equation with the linearized BGK operator that centered at bulk velocity u is: Following the classical treatment of the half-space equations, we shift the center of the Maxwellian M u to the origin by performing the change of variable v − u → v. The half-space equation (2.4) then becomes and the null space of L becomes with M is the global Maxwellian centered at the origin: It is well known that equation (2.5) is not necessarily well-posed unless appropriate boundary conditions are imposed at x = ∞. In this paper we will adopt the framework developed in [CGS88]. To illustrate the boundary conditions used in this setting, we first define a new set of basis functions for Null L. Let such that: where α, β ∈ {+, −, 0}, c = 3/2, and Using these new basis functions, we can decompose Null L into subspaces: where again β ∈ {+, −, 0}. For each fixed u ∈ R, denote the dimensions of these subspaces as Note that ν ± , ν 0 change with u. In particular, we have the following categories: As mentioned in the introduction, we will present a variational formulation using the mixed regularity idea in [ES12] to prove Theorem 2.1 for general linear or linearized half-space equations when L satisfies (A1)-(A4). We need to overcome the difficulty due to the non-coercivity of L. Although this degeneracy of L in some cases can be handled by carefully choosing appropriate function spaces for the variational formulation, we prefer to work with strictly dissipative operators.
To this end, we utilize the idea developed in [Gol08,UYY03] to modify the original equation (2.5) by adding in damping terms. The particular damping term are chosen in the way such that we can easily recover the undamped equation (2.5) for any incoming data.
The rest of this section will be devoted to the description of the variational formulation of the damped half-space equation.
2.2. Variational Formulation. We will use V = R 3 as the setting to explain the variational formulation. Other spaces for v will work in a similar way. Let u ∈ R be given. We modify (2.2) by adding in damping terms. The modified equation for f ∈ (L 2 ( dx dv)) d has the form Here the damped collision operator L d is defined by Here the constant α satisfies that 0 < α ≪ 1. More specific requirement of α will be explained in Section 5. We define the shifted "even" and "odd" parts of a function as Define the function space Thus the norm of Γ is equivalent to ( 1 a dv dx)) d . Moreover, every element g ∈ Γ has a well-defined trace: and (2.17) The inner product ·, · x,v in (2.18) is The variational formulation of (2.10) has the form Here the linear operator l(·) is given by where φ is the given incoming data and ψ + is the even (with respect to −u) part of ψ as defined in (2.13).

Main results
In this section we state the main analytical results of this work, and present the ideas of our numerical algorithms. Our strategy is to first solve the damped equation (2.10) and then recover from it the solution to the original undamped equation (2.2). Hence, in the first part we state the well-posedness of the damped equation and the quasi-optimality of the associated Galerkin method for the damped equation. In the second part we explain the details of the recovery procedures and state the convergence result for the original equation.
3.1. Solution of the damped equation.
Remark 3.1. Note that f − f ∞ for the neutron transport equations satisfy the exponential decay as x → ∞ by the comment after (2.3).
Using the variational formulation (2.19), we will solve the damped equation (3.1) by a Galerkin method.
is an orthonormal basis of (L 2 ( dv 1 )) d such that 2n (v 1 ) is even in v 1 with respect to −u for any n ≥ 1; n=1 are orthonormal bases for (L 2 ( dv 2 )) d and (L 2 ( dv 3 )) d respectively. Define the closed subspace Γ N K as Then there exists a unique f N K ∈ Γ N K such that where B and l are defined in (2.18) and (2.20) respectively. The coefficients {a kmn (x)} satisfy that The following Proposition reformulates (3.4) into an ODE with explicit boundary conditions. .
Define two 6-tensors A and B as for 1 ≤ i, k ≤ 2N + 1, 1 ≤ j, m ≤ K, and 1 ≤ l, n ≤ K. Then the variational form (3.4) is equivalent to the following ODE for the coefficients a together with the boundary conditions at x = 0: for i = 1, · · · , N and j, l = 1, 2, · · · , K.
Since our numerical examples are both in 1-dimension, we apply Proposition 3.2 to V ⊆ R 1 and obtain the following Corollary for two special cases: is an orthonormal basis of (L 2 ( dv)) d such that • ψ 2n−1 is odd and ψ 2n is even in v with respect to −u for any n ≥ 1; Define the closed subspace Γ N as Then there exists a unique f N ∈ Γ N such that where B and l are defined in (2.18) and (2.20) respectively. Moreover, the variational form (3.9) is equivalent to the following ODE for the coefficients a k (x) together with the boundary conditions at x = 0: We now state the quasi-optimality of the Galerkin approximation.
Proposition 3.5 (Quasi-Optimality). Suppose f N is an approximation defined in Proposition 3.4, and f is the unique solution to the damped equation. Then there exits a constant C 0 such that where · Γ is the norm defined in (2.14).
Numerically, the approximate solutions f N K in (3.6) (or f N (3.10) in 1D) will be solved using the method of generalized eigenvalues. In particular, we define the generalized eigenvalues and its associated eigen-tensor for (A, B) as λ ∈ R and η = (η kmn ) (2N +1)×K×K such that When reduced to 1D system, the generalized eigenvalue problem for (A, B) becomes To solve for the coefficient a(x), we take (3.10) as an example. Define γ(x) = η T B a(x) and multiply (3.10) by η T from the left. We then obtain the equation for γ as If λ = 0, then we immediately get the constraint If λ = 0, then we have γ(x) = e x/λ γ(0) .
Depending on the signs of the eigenvalues, γ either grows exponentially to infinity or decays exponentially to zero; as we look for bounded decaying solutions, this gives us constraints to γ(0) for the growing modes: If λ > 0, then we have the constraints Note that we do not need constraints for modes with negative eigenvalues. The total number of constraints in the form of (3.16) is determined by the number of positive generalized eigenvalues. The following Proposition gives the signature of (A, B): Proposition 3.6. Let A, B be the 6-tensors defined in (3.5) with any arbitrary u ∈ R and N, K ≥ 1. Then there are N K 2 positive generalized eigenvalues, N K 2 negative eigenvalues, and K 2 zero eigenvalue for the pair (A, B).
A direct application of Proposition 3.6 to the 1D case gives: Corollary 3.7. Let A, B be the matrices defined in (3.12) with any arbitrary u ∈ R and N ≥ 1. Then there are N positive generalized eigenvalues, N negative eigenvalues, and one zero eigenvalue for the pair (A, B).
Therefore in total we have N + 1 equations for a(0) given by the constraints (3.15) and (3.16). Combining them with the N equations given by the boundary conditions (3.11) for a(0), we get 2N + 1 equations for 2N + 1 unknowns {a k (0)}. The linear system (3.10) for a is then uniquely solvable, which further uniquely determines the approximate solution f N (x, v) by (3.8).

Recovery of the undamped solution.
From the solution of the damped equation (2.12), we can explicitly construct solutions to the original undamped equation (2.2). Specifically, for each 1 ≤ i ≤ ν + , let g +,i be the solution to (2.12) with boundary conditions given by X +,i : Similarly, for each 1 ≤ j ≤ ν 0 , denote g 0,j as the solution to (3.1) where the incoming boundary data is given by X 0,j . Let C be the block matrix defined by where U f = u +,1 , · · · , u +,ν+ , u 0,1 , · · · , u 0,ν0 Proposition 3.8 (Recovery). Let φ ∈ (L 2 (a(v)1 v1+u>0 dv)) d . Then C is invertible, hence (3.18) is uniquely solvable. Moreover, is the unique solution to the half-space equation where f φ,∞ ∈ H + ⊕ H 0 is the end-state given by Applying the recovery construction in the finite dimensional numerical space Γ N , we obtain the numerical solution of the undamped equation (2.2) in Γ N . The resulting approximation for the undamped equation is almost quasi-optimal with a correction term.
Proposition 3.9. Suppose f φ,N is constructed as in Proposition 3.8 with f , g +,i , g 0,j being numerical approximations obtained in Proposition 3.4 with appropriate boundary conditions. Suppose f φ is the unique solution to the equation (2.5). Then there exists a constant C 0 such that where · Γ is the norm defined in (2.14) and Remark 3.2. Note that in the above reconstruction scheme, the solutions g +,i for 1 ≤ i ≤ ν + and g 0,j for 1 ≤ j ≤ ν 0 can be precomputed, as they do not depend on the prescribed incoming data φ. In particular, we can use a higher order approximation (larger N ) for these functions.

Convergence Proof
In this section we will show the proofs for the Propositions stated in Section 3. 4.1. Inf-sup condition and well-posedness. In order to prove the well-posedness of the weak formulation (2.19), we show that the bilinear operator B satisfies an inf-sup condition. More precisely, we have Proposition 4.1 (Inf-sup). Let Γ and B be the function space and the bilinear operator defined in (2.14) and (2.18) respectively. Then B : Γ × Γ → R satisfies that Proof. Note that B is symmetric in its variables. Hence it suffices to show that the second condition in (4.1) holds. To this end, let f ∈ Γ be arbitrary. We only need to find an appropriate ψ such that Indeed, if ψ satisfies (4.2), then one can simply let Ψ = ψ ψ Γ and obtain the second inequality in (4.1) (with a different constant). The construction of such ψ will be carried out in two steps. First, let ψ 1 = f . Then Write where γ ± 's are defined as In addition, if ν 0 = 0, then (4.5) Since the matrix (v 1 + u)L −1 ((v 1 + u)X 0,k ), X 0,m v is strictly positive, there exists a constant c 0 > 0 such that Hence by multiplying (4.3) by a large enough number and adding it to (4.5), we have (4.7) B(f, ψ 1 ) ≥ κ 0,1 f 2 (L 2 (a dv dx)) d for some κ 0,1 > 0 , provided 0 < α ≪ 1. Next, let We claim that ψ 2 ∈ Γ. Indeed, by the definition of a(v), one can find two constants c 1 , c 2 > 0 such that .
Here the constants c 1 , c 2 depend on u. Thus ψ 2 ∈ (L 2 (a dv dx)) d because Moreover the definition of ψ 2 implies that Hence ψ 2 ∈ Γ and it satisfies (4.8) Using ψ 2 in B, we have for some constant κ 2 > 0. Hence by taking κ 3 > 0 large enough, we have that for some κ 0 > 0. Recall that by the definition of ψ 1 and (4.8), we also have which, together with (4.9), shows the inf-sup property of B on Γ × Γ. Next we show that if a(v) = 1 + |v|, then there exists β > 0 such that f ∈ (L 2 (e 2βx dx; L 2 (a dv))) d . The proof will be along the same line for the general case of a(v). We use the standard way to incorporate the exponential into the bilinear form by changing f by e −βx f . The new bilinear form B β is where B(f, ψ) is defined in (2.18). Note that by Cauchy-Schwartz, if we choose 0 < β ≪ α, then by the spectral gap assumption (A4), we have Hence, this extra β-term will not affect the inf-sup estimate.
Given the inf-sup condition of the bilinear form B, we can now prove the well-posedness of the weak formulation (2.19) by using the Babuška-Aziz lemma [BA72]. There are two parts in this lemma and we recall its statement below.
Moreover, f N gives a quasi-optimal approximation to the solution f in (a), that is, there exists a constant κ 1 such that Now we proceed to the proof of Proposition 3.1.
Proof of Proposition 3.1. It is straightforward to verify the boundedness of B and l defined in (2.18) and (2.20). The well-posednes of the variational form is then an immediate consequence of Proposition 4.1 and part (a) of the Babuška-Aziz lemma. In order to show that (v 1 + u)∂ x f ∈ (L 2 ( 1 a dv dx)) d , we note that the damped equation (3.1) holds in the sense of distributions by choosing the test function ψ ∈ C ∞ c ((0, ∞) × R). ∞). Therefore, if we choose φ ∈ C ∞ [0, ∞) and integrate by parts in the variational form (2.19), then boundary terms satisfy

By the density argument this implies that
(v 1 + u) φ ψ + dv at x = 0 .

4.2.
The generalized eigenvalue problem. Let us now prove the properties of the generalized eigenvalue problems (3.13) and (3.14).
Proof of Proposition 3.6 and Corollary 3.7. We first verify that (A, B) has N positive, N negative, and one zero generalized eigenvalues. By the definition of B and the strict coercivity of L d , the matrix B is symmetric and strictly positive definite. Hence the numbers of positive, negative, and zero generalized eigenvalues are the same with the signature of the matrix B −1 A. Furthermore, by the Sylvestre's Law of Inertia, B −1 A and A have the same signature. Hence, we only need to count the numbers of positive, negative, and zero eigenvalues of A. Note that by the definition of the basis functions ψ k in (5.4), A is independent of u since one can perform a change of variable v + u → v in each entry in A. Thus we only need to study the matrix A 0 with u = 0. Change the order of the basis functions such that ( ψ 1 , ψ 2 , · · · , ψ N +1 , ψ N +2 , · · · ψ 2N +1 ) = (ψ 1 , ψ 3 , · · · , ψ 2N +1 , ψ 2 , · · · , ψ 2N ) = P (ψ 1 , ψ 2 , · · · , ψ 2N +1 ) , where P is the similarity matrix. Defined A 0 = P A 0 P −1 . Then A 0 and A 0 have the same signature. By the even/odd properties of ψ i , the matrix A 0 has the form . Suppose η = (η 1,1 , · · · , η 1,N , η 2,1 , · · · , η 2,N +1 ) T = (η T 1 , η T 2 ) T is an eigenvector of A 0 with eigenvalue λ. Then one has It is clear that (η 1 , −η 2 ) is also an eigenvector of A 0 and the associated eigenvalue is −λ. This shows the eigenvalues of A 0 appear in pairs. Since A 1 has a full rank N , we have that rank A 0 = 2N . Therefore A 0 , thus A 0 and A, has N positive eigenvalues, N negative eigenvalues, and one zero eigenvalue. Now we claim that each generalized eigenpair (λ, v) of A gives rise to K 2 eigenpairs of A. Indeed, let {w (m) } K m=1 be a set of basis vectors of R K . Choose the 3-tensor η (mn) = v ⊗ w (m) ⊗ w (n) . Then for any 1 ≤ m, n ≤ K. Thus each (λ, v ⊗ w m ⊗ w (n) ) is an eigenpair of A. Note that we can also view A and B as two matrices of size ((2N + 1)K 2 ) × ((2N + 1)K 2 ) by defining a bijection between the indices Υ : {(i, j, l)| i = 1, · · · , 2N + 1, j, l = 1, · · · , K} → {1, · · · , (2N + 1)K 2 } .
Then B is symmetric and positive definite and A is symmetric. Therefore, by a similar argument as for (A, B) using Sylvestre's Law of Inertia, the number of positive, negative, and zero generalized eigenvalues agree with the signature of A. This shows there are N K 2 positive, N K 2 negative, and K 2 zero generalized eigenvalues for (A, B).

Recovering the Undamped Equation.
In this part we explain the procedures to recover the undamped equation from the damped one, thus proving Proposition 3.8.
Proof of Proposition 3.8. First we recall [CGS88] the uniqueness property of the solution to (3.19): if f is a solution to (3.19) which satisfies f ∈ (L 2 (a dv dx)) d and (v 1 + u)∂ x f ∈ (L 2 ( 1 a dv dx)) d , then f must be unique. For the convenience of the reader, we brief explain its proof: Suppose h is a solution to the half-space equation (3.19) with incoming data φ = 0. Then R 3 (v 1 + u)h 2 dv is decreasing in x. Since there exists h ∞ ∈ H + ⊕ H 0 such that h − h ∞ ∈ (L 2 ( dv dx)) d , we can find a sequence x k such that This holds in particular at x = 0. Since the incoming data is zero at x = 0, the outgoing data at x = 0 must also be zero and The conservation property of the half-space equation then implies that h(x, ·) ∈ (Null L) ⊥ . By multiplying the equation by h and integrate over v, we have h, Lh = 0 for all x ≥ 0. Hence Ph = 0 for all x ≥ 0 by the spectral gap of L in (A4). Therefore h ≡ 0 and the solution to the half-space equation is unique. Denote and (4.10) We will separate two different cases according to dim H 0 .
Case 1: dim(H 0 ) = 0. In this case, there are only various numbers of X + 's and X − 's. Hence (3.1) reduces to (4.11) and U reduces to Multiplying (4.11) by X +,k , X −,j and integrating over v ∈ R, we obtain a linear system for U : where the coefficient matrix is diagonal: where D + , D − are positive definite if they exist and D + = diag(γ +,1 , · · · , γ +,ν+ ) , where γ ±,k are defined as in (4.4). Since solutions to (4.11) are in (L 2 ( dv dx)) d , it is clear that Hence (4.11) further reduces to (4.14) Given any φ ∈ L 2 , suppose f is the solution to (4.14) (or (3.1)) with φ as the incoming data. Let g k be the solution to (4.14) (or (3.1))with the incoming data We claim that (4.15) (v 1 + u)X +,k , g k v x=0 = 0 , for each k ≥ 1 .
Note that D 2 is symmetric. Hence, D 2 has a complete set of eigenvectors. Let (λ, E) = (λ, (e 1 , e 2 ) T ) be an eigenpair of A 4 . Then which is equivalent to e 2 = λe 1 , e 1 T + αD 2 e T 2 = λe T 2 . Note that λ = 0. Since D 2 is symmetric, it has a complete set of orthogonal eigenvectors. Let e be an arbitrary eigenvector of D with eigenvalue λ e and take e 2 = e. Then (4.23) e 1 = 1 λ e , 1 λ e T + αD 2 e T = λe T .

Thus
(4.24) 1 λ + αλ e − λ = 0 , which has exactly one negative solution for λ. Since the set of eigenvectors of D 2 is complete, the matrix A 4 has exactly ν 0 negative eigenvalues. Together with D − , we have that A 3 , thus A 2 , has exactly ν − + ν 0 negative eigenvalues, which prescribes ν − + ν 0 conditions on U such that where E k are the eigenvectors associated with negative eigenvalues. Write each E k as E k = (e k,+ , e k,− , e k,0 , e k,L,0 ) T . Now we make the following claim: Let QN 2 as the space spanned by the eigenvectors of A 2 with negative eigenvalues. Then there exists a nontrivial F ∈ QN 2 such that Since QN 2 is an invariant subspace of A 2 , we have that where f 2 = αA T 21 f T + + (I + αB) 1/2 (αB) 1/2 ) T f T 0 . By the symmetry and non-degeneracy of A 2 , the quadratic form given by A 2 on QN 2 is strictly negative. Therefore, Since D + is strictly positive definite, we have that f + = 0 and which implies that F = 0. This contradicts the assumption that F is non-trivial. Hence the matrix E is non-singular.
Now we are ready to show the recovery of solutions to the undamped equation from (3.1). Let g +,j , g 0,k be the solutions to (3.1) with the incoming data X +,j , X 0,k with 1 ≤ j ≤ ν + and 1 ≤ k ≤ ν 0 respectively. Let C be the matrix defined in (3.17).

4.4.
Quasi-Optimality. We now prove Proposition 3.2, 3.5, and 3.9. The key observation in [ES12] for the neutron transport equation is that the inf-sup condition holds on the finite-dimensional space as long as (v 1 + u)∂ x f + stays in this subspace. In this case, the proof for the well-posedness of the finite-dimension system will be almost a verbatim as for the original equation. Here we have a similar situation for the damped equation (3.1) with a slight modification.
Proof of Proposition 3.2 and 3.5. For the well-posedness part, we only explain the modification in choosing the test functions ψ 1 and ψ 2 . For any f ∈ Γ N , we choose where P N : (L 2 ( dv)) d → Γ N is the projection onto Γ N . The rest of the estimates are similar to the proof in Section 5, thus omitted. Now we show that the boundary conditions for the solution to (3.9) are given by (3.11). To this end, we first choose test functions G 2j (x, v) = g(x)ψ 2j (v) where g(x) ∈ C ∞ c ([0, ∞)) and 1 ≤ j ≤ N . Applying G j in (3.9), we get By integration by parts in (4.32) we obtain for each 1 ≤ j ≤ N . The set of N boundary conditions (3.11) then follows from integrating by parts in (4.34) and applying (4.33).
The error estimate of the undamped equation follows from Proposition 3.5.
Proof of Proposition 3.9. First we note that there exist constants κ 4 , κ 4 > 0 such that where Φ is defined in (4.18) or (4.28), Φ N is given by the same formula for Φ with f changed to f N , and g +,i,N , g 0,j,N are the approximations in Γ N to g +,i , g 0,j . Second, since f − g is a solution to the damped equation, we have where f N , g N ∈ Γ N . Therefore, Note that the second inequality holds because Φ ∈ H + ⊕ H 0 ⊆ Γ N .

Numerical Examples
As explained in Section 3, the overall strategy to solve the half-space equation consists of two steps: First, we solve for numerical solution to the half-space damped equation (3.1) using the Galerkin approximation; Second, we recover the undamped solution by Proposition 3.8, which involves the solution of the damped equation with various boundary conditions to obtain the matrix C in the linear system (3.18).
We will consider the linearized BGK equation and the linear transport equation below. We will restrict ourselves here to one dimensional examples. As in Proposition 3.4, for the Galerkin approximation, we specify a set of even and odd functions to form the approximation space Γ N , whose choice depends on the particular equation under study. Using Corollary 3.7, the solution of the system under Galerkin approximation (3.10)-(3.11) is reduced to solving the generalized eigenvalue problem (3.14), for which we need to assemble the matrices A and B using Gaussian quadrature. This will be discussed in more details below.
Our algorithm is implemented in MATLAB. The Gaussian quadrature abscissas and weights are obtained using symbolic calculations to guarantee precision. The orthogonal polynomials can be constructed using three term recursion formula (see for example [Shi81]), and the details are recalled in Appendix A for completeness.
The basis functions ψ we need are either odd or even with respect to v = −u, hence we shift the functions B n 's by −u and make even / odd extensions: Finally, the basis functions ψ k 's are obtained by multiplying these functions by the square root of the Maxwellian: for n ≥ 1 By definition, ψ 2n−1 is odd, ψ 2n is even, and they form a orthonormal basis of L 2 (dv). For a fixed n, (v + u)ψ 2n (v) is a odd function with respect to v = −u. For v > −u, Since (v + u)B n (v + u) is a n-th order polynomial in v + u, there exists an expansion This yields that Therefore, Γ N = span{ψ 1 , · · · , ψ 2N +1 } satisfies the condition of Proposition 3.4 and the variational formulation (3.9)-(3.11) is well-posed. The (2N + 1) × (2N + 1) matrices A and B are then given by Note that both matrices are symmetric. The matrix A can be obtained by using the recurrence relation of the orthogonal polynomials. For the matrix B, recall that All the integrals involved in calculating B can be easily made exact up to machine precision by using Gaussian quadrature. For simplicity, let us just focus on ψ 2j χ 0 dv and note that the other integrals share the same structure, that is, the integrand is a product of two polynomials and two Gaussians e −v 2 /2 and e −(v+u) 2 /2 . To evaluate this type of integral using Gaussian quadrature, we first split the integral into two parts: Note that ψ 2j , on either side of −u, is a (j − 1)-th order polynomial multiplied by exp(−(v + u) 2 /2), while χ 0 is a quadratic function multiplied with a different weight function exp(−v 2 /2). The product of two Gaussians centered at different locations could be combined into a single Gaussian: Similarly, we have the integration for v < −u: The integrals (5.7) and (5.8) can be evaluated up to machine precision by Gaussian quadrature based on weight e −(v−u/2) 2 and e −(v+u/2) 2 respectively, as B j−1 χ 0 e v 2 /2 is a polynomial with its degree up to N + 3. The boundary condition (3.11) requires the numerical evaluation of the integral v+u>0 (v + u)φψ 2j dv.
We calculate this using Gaussian quadrature with the weight e −(v+u) 2 . The error of the quadrature depends on the number of quadrature points and the regularity of the incoming data φ.
We now present some numerical results for the linearized BGK equation. In the first set of examples, we compare our numerical results with analytical solutions, when the specified boundary data φ is given by the restriction of some f ∈ H 0 ⊕ H + on v > −u. In this case, the solution to the undamped equation (1.1) is simply f on the whole velocity space. As discussed in (2.9), the dimension of the space H 0 ⊕ H + depends on the bulk velocity u and the sound speed, which is c = 3/2 in our case as T = 1/2. We will choose χ +/−/0 defined in (2.7) as the incoming data. By the uniqueness of the half-space equation, the solution will simply be χ +/0 when the incoming data is chosen as χ +/0 . We take six choices of u corresponding to the six cases listed in (2.9) (the case u < −c gives an empty H 0 ⊕ H + hence not included). The results are shown in Figures 1-6 below. In all these figures, the blue squared line is the incoming data, given by χ − , χ 0 and χ + respectively. The green triangle line is the solution at x = ∞, and the red dotted line is the solution at x = 0.
Several remarks are in order: First, when the χ modes lie in H 0 ⊕ H + for the given bulk background velocity u, we observe in Figure 1-6 that the solution at x = 0 gives a perfect match. We thus recover the exact solution from the numerical scheme. Second, we note that in general, the solution exhibits a jump at v = −u, as clearly seen for example in Figure 1(left). This justifies our choice of the even-odd formulation and basis functions from the half-space Hermite polynomials. Finally, we remark that we have used a filtering (with 2nd order cosine filter) to reduce the Gibbs oscillations caused by the large derivatives in some cases (for instance Figure 2(left)).  Next, we consider an example where the exact solution is not known. We solve the equation (1.1) for u = 0 with boundary data φ = v 3 , v > 0. The numerical solution is shown in Figure 7.
5.2. Isotropic neutron transport equation. We further consider the isotropic neutron transport equation. The construction of the basis functions is similar to the linearized BGK case. However, instead of using half-space Hermite polynomials, we start with Legendre polynomials on the interval [0, 1] and carry out even / odd extensions. The Legendre polynomials, which are orthogonal polynomials for constant weight   function, are used as the equilibrium states for the neutron transport equation are simply constants. We then apply Gauss-Legendre quadrature to assemble A and B for the generalized eigenvalue problem. We will skip further details here, as the construction is relatively straightforward compared with the linearized BGK case.
To validate our methods in this case, we compare numerical solution with analytical solution with boundary data given by φ = v for v ∈ [0, 1]. The analytical solution is known in case as    where H is the Chandrasekhar H-function. In Figure 8 we plot both analytical and numerical solutions. The plot shows good agreement of the numerical solution with the exact one. Furthermore, the limit at x = ∞ of the solution to the half-space isotropic NTE is a constant, whose amplitude is known as the extrapolation length. In Table 1 we compare our numerical approximation of extrapolation length with the exact result, which is again in good agreement. In comparison, we note that the approximate value for the extrapolation length obtained in [Cor90] is 0.71040377 with 70 modes, while we achieve better results with piecewise polynomial of only 12-th order.  Integrate the identity (A.5) over v with the weight, we get where the second equality is obtained by taking the inner product with B ′ n+1 of recursion equation (A.2), and the third comes from integration by parts. From this we get the first deduction relation in (A.3). Next multiply (A.5) with v and then integrate, we obtain where the first equality comes from the fact that ∞ 0 vB n+1 B ′ n e −(v−u) 2 dv = 0, the second is due to integration by parts, and the third comes from integrating the recursion equation (A.2) multiplied by vB n+1 . This gives the other deduction relation in (A.3).