Weak Solutions of the Relativistic Vlasov-Maxwell System with External Currents

The time evolution of a collisionless plasma is modeled by the relativistic Vlasov-Maxwell system which couples the Vlasov equation (the transport equation) with the Maxwell equations of electrodynamics. We consider the case that the plasma consists of $N$ particle species, the particles are located in a bounded container $\Omega\subset\mathbb R^3$, and are subject to boundary conditions on $\partial\Omega$. Furthermore, there are external currents, typically in the exterior of the container, that may serve as a control of the plasma if adjusted suitably. We do not impose perfect conductor boundary conditions for the electromagnetic fields, but consider the fields as functions on whole space $\mathbb R^3$ and model objects, that are placed in space, via given matrix-valued functions $\varepsilon$ (the permittivity) and $\mu$ (the permeability). A weak solution concept is introduced and existence of global in time solutions is proved, as well as the redundancy of the divergence part of the Maxwell equations in this weak solution concept.


Introduction
boundary conditions on ∂Ω have to be imposed. In the exterior of Ω, there are external currents, for example in electric coils, that may serve as a control of the plasma if adjusted suitably. In order to model materials that are placed somewhere in space, for example (almost perfect) superconductors, we consider the permittivity ε and permeability µ, which are functions of the space coordinate, take values in the set of symmetric, positive definite matrices of dimension three, and do not depend on time, as given. With this assumption we can model linear, possibly anisotropic materials that stay fixed in time. We should mention that in reality ε and µ will on the one hand additionally depend on the particle density inside Ω and on the other hand additionally locally on the electromagnetic fields, typically via their frequencies (maybe even nonlocally because of hysteresis). However, this would cause further nonlinearities which we avoid in this work.
The unknowns are on the one hand the particle densities f α f α (t, x, v), α 1, . . . , N, which are functions of time t ∈ R, the space coordinate x ∈ Ω, and the momentum coordinate v ∈ R 3 . Roughly speaking, f α (t, x, v) indicates how many particles of the α-th species are at time t at position x with momentum v. On the other hand there are the electromagnetic fields E E(t, x), H H(t, x), which depend on time t and space coordinate x ∈ R 3 . The Dand B-fields are computed from E and H by D εE and B µH. We will only view E and H as unknowns in the following. The main assumption about ε (and likewise µ) in Section 3 will be σ ≤ ε ≤ σ ′ for some σ, σ ′ > 0 in the sense of positive definiteness. This property implies that is a norm on L 2 R 3 ; R 3 , which is equivalent to the standard L 2 -norm.
The Vlasov-Maxwell system on a time interval with given final time 0 < T • ≤ ∞, equipped with boundary conditions on ∂Ω and initial conditions for t 0, is then given by the following set of equations; we explain the appearing notation afterwards:  H) at time t 0, that is to say the function f α (0, ·, ·). We will use this notation often, also similarly for other functions. Note that throughout this work we use modified Gaussian units such that the speed of light is normalized to unity and all rest masses m α of a particle of the respective species are at least 1. In (VM.1), e α is the charge of the α-th particle species and v α the velocity, which is computed from the momentum v by v α v Clearly, v α < 1, that is, the velocities are bounded by the speed of light. Equation (VM.2) describes the boundary condition on ∂Ω. Typically, one imposes specular boundary conditions. Thus it is natural to consider the following decompositions: where n(x) is the outer unit normal of ∂Ω at x ∈ ∂Ω and 0 < T ≤ ∞. In (VM.2), f α ± are the restrictions of f α to γ ± T • . The operator K α maps functions on γ + T • to functions on γ − T • . In Section 3 we deal with the case that describes reflection on the boundary and a α , satisfying 0 ≤ a α ≤ 1, describes how many of the particles hitting the boundary at time t at x with momentum v are reflected (and not absorbed); g α is the source term according to how many particles are added from outside. We will deal with partially absorbing (a α ≤ a α 0 for some a α 0 < 1) and purely reflecting (a α 1, g α 0) boundary conditions. In (VM.4) and (1.2a), j and ρ are the current and charge density. Typically they are the sum of the internal current and charge densities, and some external current density u and charge density ρ u resulting from u. Usually, the divergence equations (1.2) are known to be redundant if all functions are smooth enough, local conservation of charge is satisfied, i.e. ∂ t ρ + div x j 0, . Some notation and (1.2) holds initially, which we then view as a constraint on the initial data. Therefore, in the first sections we ignore (1.2) and discuss in Section 4 in what sense (1.2) is satisfied in the context of a weak solution concept. Since (1.2) has to hold on whole space R 3 , the main difficulty will be that we have to "cross over" ∂Ω.
The paper is organized as follows: In Section 2.3 we state our main two theorems. The first regards the existence of weak solutions to (VM). In Section 3 we prove this theorem. To this end, we state some basic results about linear Vlasov and Maxwell equations (Section 3.1), approximate the given functions in a proper way (Section 3.2), consider a cut-off system (Section 3.3), and finally remove the cut-off (Section 3.4). The second main result regards the redundancy of the divergence equations in our weak solution concept. We prove this theorem in Section 4 and give some comments on the physical interpretation of the obtained equations.
In the first part, we proceed similarly to Guo [11], who proved existence of weak solutions in the case that ε µ 1, u 0, and the electromagnetic fields are subject to perfect conductor boundary conditions on ∂Ω, i.e., E × n 0. However, there is no need of artificially inserting the factor e −t as is done throughout that paper. The more important motivation of our paper is the following: The papers concerning plasma in a domain we are aware of deal with perfect conductor boundary conditions for the electromagnetic fields. Such a set-up can model no interaction between this domain and the exterior. However, considering fusion reactors, there are external currents in the exterior, for example in field coils. These external currents induce electromagnetic fields and thus influence the behavior of the internal plasma. Even more important, the main aim of fusion plasma research is to adjust these external currents "suitably". Thus, we impose Maxwell's equations globally in space and model objects like the reactor wall, electric coils, and almost perfect superconductors via ε and µ.
Vlasov-Maxwell systems have been studied extensively. In case of no reactor wall, i.e., the Vlasov equation is imposed globally in space (as well as Maxwell's equations), global well-posedness of the Cauchy problem is a famous open problem. Global existence and uniqueness of classical solutions has been proved in lower dimensional settings, see Glassey and Schaeffer [5,6,7,8]. In the full three-dimensional setting, global existence of weak solutions was proved by Di Perna and Lions [2]. Their momentum-averaging lemma is fundamental for proving existence of weak solutions in any setting (with or without boundary, with or without perfect conductor boundary conditions and so on), since it handles the nonlinearity in the Vlasov equation. However, uniqueness of these weak solutions is not known. For a more detailed overview we refer to Rein [13].

Some notation
In the following, we denote by χ M the characteristic function of some set M (i.e., χ M (x) 1 if x ∈ M and 0 otherwise) and by χ T the characteristic function of [0, T]. For 1 ≤ p < ∞ and α 1, . . . , N we define . Weak formulation equipped with the corresponding weighted norm. Here, A ⊂ R 3 × R 3 or A ⊂ R × R 3 × R 3 is some Borel set equipped with a measure a and the weight v 0 α is given by where 0 ∈ I ⊂ [0, ∞[ is some interval, G is some C k or L p , and X is a normed vector space.
For ease of notation it will be convenient to introduce a surface measure on [0, ∞[×∂Ω×R 3 , namely Since ε is already used for the permittivity, the letter ι, and not ε, will always denote a small positive number.
For a matrix A ∈ R 3×3 and a positive number Finally, for a normed space X, some x ∈ X and r > 0, B r (x) denotes the open ball in X with center x and radius r. Furthermore we abbreviate B r : B r (0).

Weak formulation
The space of test functions for (VM.1) to (VM.3) will be On the other hand, will be the space of test functions for (VM.4) to (VM.6). We start with the definition of what we call solutions to (VM).

. Statement of main results
We call a tuple f α , f α + α , E, H, j a weak solution of (VM) on the time interval I T • if (for all α) We easily derive this weak formulation after multiplying the respective equations of (VM) with the respective test function and integrating by parts, assuming all functions are smooth enough.

Statement of main results
We have two main results: The first is about existence of weak solutions in the case of partially absorbing boundary conditions for particle species 1, . . . , N ′ and purely reflecting boundary conditions for particle species N ′ + 1, . . . , N. We assume that the following conditions hold:

Statement of main results
Estimate on j int : (2.8) The second main result answers the question whether the divergence equations (1.2) are automatically satisfied if we have a weak solution of (VM). To this end, we have to introduce an external charge density ρ u corresponding to u and assume that local conservation of the external charge holds: which is to be understood in the following weak sense: Here, ρ u andρ u are extended by zero outside Γ.
Then our second main result is (see Section 4): such that the tuple f α , f α + α , E, H, j int + u is a weak solution of (VM) in the sense of Definition 2.1. Furthermore, assume that Condition 2.4 holds. Moreover, let initially div x εE 4π ρ int +ρ u : 4π (iii) If, additionally to the given assumptions, f α is given by Note that K α need not necessarily have the structure (1.3) in Theorem 2.5.

Existence of weak solutions
In this section, we proceed similarly to Guo [11] with necessary modifications being made, who considered the problem with ε µ 1, u 0, and perfect conductor boundary conditions for the electromagnetic fields on ∂Ω. Citations of this paper always refer to the relativistic version of the respective lemma, theorem etc., see [11,Section 5].

Results about linear Vlasov and Maxwell equations
The strategy is to consider an iteration scheme where we decouple Vlasov's equations from Maxwell's equations in each iteration step and hence only have to solve linear problems. Thus it is natural to consider linear Vlasov and Maxwell equations first. Regarding the Vlasov part, we refer to Beals and Protopopescu [1]. Considering the linear problem (on with a Lipschitz continuous, bounded force field F, that is divergence free with respect to v, they introduced a space of test functions associated to F. As in [11, Lemma 2.1.] we can show that our test function space Ψ T belongs to that test function space for each F and T, where one needs the assumption that ∂Ω be of class C 1,κ and that the support of any ψ ∈ Ψ T be away from γ 0 T and {0} × ∂Ω × R 3 . In [1], "strong" solutions in a set of L p -functions for which a trace on the boundary exists in the sense of the following extended Green's identity were searched for: which is supposed to hold for all test functions φ. Here, D ± T are the outgoing/incoming sets associated to the characteristic flow of Y and dν ± are associated measures. In our case, we can split D + [1]). Then, dν ± dγ α on γ ± T and dν ± dvdx on t 0 and t T, and we decompose Lipschitz continuous, bounded, and divergence free with respect to v, and letf ∈ L 1 ∩ L ∞ Ω × R 3 ,

. Results about linear Vlasov and Maxwell equations
for any 0 < T ∈ I T • and 0 < R < ∞.
Proof. By [1, Theorem 1], there is a unique, strong solution of (3.1) for each T ∈ I T • . Since using the convexity of the p-th power. This yields , dγ α and the respective norms coincide.
To prove the second estimate, let Noticing that Y β f F · β ′ f and using the 1-norm balance of [1, Proposition 1] we get by

Results about linear Vlasov and Maxwell equations
Writing the terms explicitly and using the fact that v 0 α is monotonically increasing in |v|, we arrive at (3.3).
where we optimize r : Regarding the linear Maxwell part on I T • , there holds the following basic result: and Proof. For the existence theory (and a definition of uniform local Sobolev spaces H k ul ) we refer to [12]. Equation (3.7) is derived straightforwardly by differentiating both sides and using the symmetry of ε and µ. We then get (3.8) by applying Lemma 3.3 using the uniform positive definiteness of ε and µ.
Here and later, we need the following version of the quadratic Gronwall lemma, which is a slight improvement of [3, Theorem 5]: Assume that the following inequality holds for all t ∈ [a, b]: Then we have Proof. Let ι > 0 and consider By assumption we have Integrating this estimate from a to t yields Since ι > 0 is arbitrary, the proof is finished.

Approximations of the data
Throughout this section we assume that Condition 2.2 is satisfied. We have to modify the data as follows to be able to apply the statements of Section 3.1: For α 1, . . . , N we define a α k : a α and for α N ′ + 1, . . . , N we define a α k : k k+1 a α . Hence all a α k are bounded away from 1. Furthermore, choose approximating sequences Additionally, we have to smooth ε and µ. In the following, have in mind that for a symmetric, positive definite matrix A ∈ R 3×3 and some Approximations of the data where we use the norm where the last equality holds for symmetric, positive definite A. Thus, for some measurable We want to construct sequences of smooth ε k , µ k with σ ≤ ε k , µ k ≤ σ ′ in such a way that these sequences converge to ε, µ in a certain sense. We perform the construction of (ε k ), the one for µ k works totally analogously.
Finally define ε k : ω s k * ε k + σI 3 . Note that ε k is smooth and constant for |x| large (and hence of class H 3 Furthermore, for any E, x ∈ R 3 it holds that Note that for the last line we used the fact that the integral of ω s over whole R 3 equals 1 for any s > 0.
. A cut-off problem

A cut-off problem
In order to construct a weak solution of (VM), we first turn to a cut-off problem where we consider bounded time and momentum domains. Whereas the cut-off in time is no real drawback, the cut-off in momentum space is on the one hand unpleasant, but on the other hand necessary. To understand this necessity, we should recall (3.8). Consider there j to be the sum of some external current and the current j int induced by the particle densities. In an iteration scheme we would like to have an estimate like (3.8) for the fields where the right hand side is uniformly bounded along the iteration. Then we could extract some weakly converging subsequence. However, for this uniformity, we would need that j int is uniformly bounded in L 1 0, T; L 2 R 3 ; R 3 along the iteration. This would require a better estimate than (3.4) where we only can put our hands on the L 4 3 R 3 ; R 3 -norm of j int (at each time).
Moreover, in an energy balance along the iteration, the crucial terms describing the energy transfer due to the internal system will not cancel out; this would only be the case if we solve (VM) simultaneously along an iteration. Now if we consider a cut-off problem (the cut-off referring to momentum space) we can simply estimate the L 1 -norm of j int by the L 2 -norm in momentum space and then use (3.2) for p 2, so we get the desired uniform boundedness along the iteration. Later, adding the limit versions of (3.3) and (3.7), we observe that the problematic terms on the right hand side, that is to say the terms ±E · j int , cancel out. Thus, now (after a Gronwall argument) having a full energy estimate with only expressions of the data on the right hand side, we find that a posteriori the cut-off does not substantially enter this estimate, so we will be able to get a solution of the system without a cut-off by considering a sequence of solutions due to larger and larger cut-off domains.
We differ from [11] as follows: Firstly, we do not have to cut off Ω, since we only consider a bounded Ω. Secondly, we solve the linear Vlasov equation on whole momentum space R 3 and not only on a cut-off domain. Our cut-off only appears in the definition of the internal current j int k . Thirdly, as already said in the introduction, there is no need of the factor e −t , and without this factor the estimates are more "natural".
To make things more precise, let 0 < R < ∞, define R * : min{R, T • }, and start the iteration . We assume that we already have iterates of the k-th satisfying E k , H k ∈ L ∞ 0, R * ; L 2 R 3 ; R 3 ∩ C 0,1 [0, R * ] × Ω; R 3 . We first solve the Vlasov part with given force field F α k : e α E k + v α × µ k H k , which is Lipschitz continuous and bounded on [0, R * ]×Ω×R 3 , and divergence free with respect to v. Indeed, we can solve (3.10) applying Proposition 3.1 and noticing that a α k+1 is bounded away from 1 on γ − R * . Therefore we have Next we want to solve the Maxwell part. Now the cut-off appears: We define the current where we integrate only over the cut-off domain B R rather than over the whole momentum space. Note that j int k+1 (u) is defined to be 0 outside Ω (Γ). By and f α k+1 ∈ L ∞ 0, R * ; L 2 Ω × R 3 we have j k+1 ∈ L 1 0, R * ; L 2 R 3 . In order to apply Propo- With this smoothed current as the source term in the Maxwell system we solve Indeed, applying Proposition 3.2, we see that there is a unique solution (E k+1 , H k+1 ) ∈ C 0, R * ; H 3 R 3 ; R 6 ∩ C 1 0, R * ; H 2 R 3 ; R 6 . By Sobolev's embedding theorems it holds that E k+1 , H k+1 ∈ C 0,1 [0, R * ] × Ω; R 3 . Altogether, the induction hypothesis is satisfied so that we can proceed with the next iteration step. In order to extract some weakly converging subsequence, we have to establish suitable estimates. To this end, consider (3.2) and (3.8) applied to (3.10) and (3.14): and .

A cut-off problem
Note that we need ε k (x), µ k (x) ≥ σ uniformly in x and k to get (3.16). For α 1, . . . , N ′ , (3.15) reduces to and to As for the boundary values, we have to distinct absorbing and reflecting boundary conditions. For α 1, . . . , N ′ , (3.17) yields the boundedness of f α k ,+ in any L p γ + R * , dγ α , 1 ≤ p ≤ ∞, so we may extract a subsequence that converges weakly in L p γ + R * , dγ α for 1 < p < ∞ and weakly-* in L ∞ γ + R * , dγ α to some nonnegative f α R,+ . For α N ′ + 1, . . . , N, (3.18) delivers a uniform estimate only for p ∞ so here we may extract a subsequence that only converges weakly-* to some nonnegative f α R, and for α N ′ + 1, . . . , N additionally (3.21) Next we turn to an estimate on the electromagnetic fields. To examine (3.16) further, we insert the properties of j k+1 on the right hand side to get for 0 < T ≤ R * using (3.12). The right hand side is bounded uniformly in k. Moreover, the first term on the right hand side of (3.16) is bounded uniformly in k by ε k , µ k ≤ σ ′ and the L 2 -convergence of the approximating initial data. Thus, we may extract a subsequence We now show that f α R , f α R,+ α , E R , H R , j R is a weak solution of (VM) on the time interval [0, R * ] in the sense of Definition 2.1. Clearly, all functions are of class L 1 loc . The main task is to show that we may pass to the limit in (2.1) and (2.2) applied to the iterates: We have for all ψ ∈ Ψ R * , ϑ ∈ Θ R * , and k ≥ 1 We can pass to the limit in (3.23) and (3.24): Whereas the terms including the curl are easy to handle by weak convergence of E k , H k , we have to take more care about the terms including ε k , µ k , and j k . For the first ones, let K ∈ N such that ϑ vanishes for |x| ≥ K so that we in fact only integrate over B K . For k ≥ K we have ε − ε k L 2 (BK;R 3×3 ) ≤ ε − ε k L 2 (Bk;R 3×3 ) < 1 k by (3.9) so that ε k → ε in L 2 B K ; R 3×3 . This is enough for passing to the limit in the terms including ε k since we additionally have E k ⇀ E R in L 2 [0, R * ] × R 3 ; R 3 , even strong convergence of the approximating initial data, and the boundedness of the time interval [0, R * ]. Similarly, we argue for the terms with µ k . So there only remains the term including j k . To tackle this one, we estimate where the first term on the right hand side converges to 0 for k → ∞ by construction of j k and each summand of the second term by weak convergence of the f α k . Note that for the latter limit our cut-off plays an important role since v α · ϑχ {|v|≤R} ∈ L 2 [0, R * ] × R 3 × R 3 .

. A cut-off problem
Passing to the limit in (3.22) is more complicated, especially because of the nonlinear product term including E k , H k , and f α k . The other terms are easy to handle due to weak convergence of f α k and weak (or weakly-*) convergence of f α k ,+ . The nonlinear term is handled as in [11, Proof of Lemma 3.1.] by a highly nontrivial tool, namely the momentum-averaging lemma (see [2], or [13] for a shortened proof). For this, it is important that the sequences for k ≥ 1 and any T ∈ ]0, R * ]. We consider the right hand sides of (3.25) and (3.26) further. The term including the initial data of the electromagnetic fields is bounded uniformly in k due to

. A cut-off problem
After approximating e α · α in L 2 B R ; R 3 by C ∞ c B R ; R 3 -functions and using the momentum averaging lemma again we have, up to a subsequence, Similarly, Unfortunately, this is not enough since we in fact have to consider E k−1 · j int k − E k · j k . To get hands on this term, choose ϕ 1 and choose u k ∈ C ∞ c ]0, R * [ × Γ; R 3 such that u − u k L 1 (0,R * ;L 2 (Γ;R 3 )) < 1 k .
Using these approximations and (3.11) and (3.13) we estimate where C > 0 does not depend on k since we already have a uniform bound on the E k in L ∞ 0, R * ; L 2 R 3 ; R 3 . Furthermore, h k is continuous with respect to T and by (3.28) and (3.29). Moreover, we have where C > 0 does not depend on k (and T) by the uniform boundedness of the E k in L ∞ 0, R * ; L 2 R 3 ; R 3 and (3.12) (combined with (3.17) and (3.18), respectively). Therefore Then there also holds . A cut-off problem By E k , H k ∈ C 0, R * ; L 2 R 3 ; R 3 , u k ∈ C 0, R * ; L 2 Γ; R 3 , and by differentiability of l k we can apply Lemma 3.3 and thus obtain altogether. For k → ∞, let A ⊂ [0, T ′ ] be measurable and integrate (3.33) over A. As for by weak convergence and we have a pointwise bound uniformly in T and k by (3.33). Additionally exploiting weak convergence and weak lower semi-continuity, respectively, the strong convergence of the initial electromagnetic fields, and (3.32) we may pass to the limit and conclude, since A was arbitrary, that for all T ∈ ]0, R * ], after taking T T ′ . This is exactly the energy estimate we wanted to derive since R does no longer appear on the right hand side.
Lastly, we show that, up to a subsequence, j int k ⇀ j int R in L . Removing the cut-off 3 4 for 0 ≤ T ≤ R * and the right hand side is bounded in L 4 3 ([0, R * ]) uniformly in k by virtue of (3.34). Therefore we may assume that j int k converges weakly in L 4 3 [0, R * ] × Ω; R 3 . It is easy to see that the weak limit has to be j int R . As for the desired bound, we proceed similarly to (3.4) and (3.5), respectively, sum over α, apply a Hölder estimate for the sum, and use the known estimates to get for any 0 < T ≤ R * .

Removing the cut-off
Finally we obtain a solution of (VM) on the time Interval I T • by letting R → ∞. To this end, it is crucial that the right hand sides of the obtained estimates of the previous section do not depend on R; see (3.19) to (3.21), (3.34), and (3.35). Take the sequence (R m ) m (m) m , then we see by a diagonal sequence argument that, for certain We may pass to the limit in the respective estimates to obtain (2.3) to (2.8). Passage to the limit in the weak formulation of (VM) works in the same way as in [11,Theorem 4.1.]. That the weak limit of the j int m is indeed the current density j int induced by the f α is proved in the same way as in [13,Proposition 4] exploiting the energy estimate.

The redundant divergence equations and the charge balance
In this section, we want to deduce in what sense the divergence equations (1.2) hold for a solution of (VM) in the sense of Definition 2.1. This is much more difficult than in [11,Lemma 4.2.] since we consider these divergence equations on whole R 3 instead of Ω. The weak formulation of (1.2) is Obviously, (4.1) is equivalent to (1.2) be satisfied on I T • × R 3 in the sense of distributions. For (1.2) should propagate in time, we have to demand that (1.2) holds initially as a constraint on the initial data, that is to say div εE 4πρ, div µH 0 on R 3 in the sense of distributions, or, equivalently, Now let f α , f α + α , E, H, j be a weak solution of (VM) on the time interval I T • . It is easy to see that (4.1b) holds: Define Clearly, ϑ ∈ Θ T • . Hence (2.2b) and ξ The redundant divergence equations and the charge balance and we are done. As for (4.1a), we have to exploit local conservation of charge. Consequently, we have to determine what ρ is and have to use the Vlasov equations (their weak form, more precisely). Therefore, we have to make use of (2.1) in order to put the internal charge density into play. However, the test functions there have to satisfy ψ ∈ Ψ T • but a test function of (4.1a) does not depend on v. Consequently, we, on the one hand, have to consider a cut-off in momentum space, and, on the other hand, have to show that (2.1) also holds if the support of ψ is not away from γ 0 For the latter one, the following technical lemma is useful. There and throughout the rest of this section, we assume that Ω ⊂ R 3 is a bounded domain such that ∂Ω is of class C 1 ∩ W 2,∞ . Here, ∂Ω being of class C 1 ∩ W 2,∞ means that it is of class C 1 and all local flattenings are locally of class W 2,∞ .
Then there is a sequence ψ k ⊂ Ψ T • such that for k → ∞ and there is 0 < r < ∞ such that ψ and all ψ k vanish for t ≥ r. Here, is compact (which can be achieved since the hyperplane where t 0 is smooth). By assumption about ∂Ω, for each For any x ∈ ∂Ω we choose an open set U x ⊂ R 3 such that x ∈ U x and U x ⊂⊂Ũ x (here, A ⊂⊂ B is shorthand for A bounded and A ⊂ B). Then ∂Ω ⊂ x∈∂Ω U x , whence there are a finite number of points, say x i ∈ ∂Ω, i 1, . . . m, such that ∂Ω ⊂ m i 1 U i , since ∂Ω is compact. Here and in the following, we write U i : U x i , U i : Ũ x i , and F i : F x i . Since it holds that Ω \ m i 0 U i . Now let ζ i , i 0, . . . , m, be a partition of unity on M subordinate to U i , i 0, . . . , m, i.e., the ζ i are of class C ∞ , 0 ≤ ζ i ≤ 1, supp ζ i ⊂ U i , and m i 0 ζ i 1 on M (and hence on Ω, in particular). Furthermore, let η ∈ C ∞ (R) such that 0 ≤ η ≤ 1, η y 0 for y ≤ 1 2 , and η y 1 for y ≥ 1.

The redundant divergence equations and the charge balance
Note that the rows are orthogonal and have length one, and that A i is of class C ∩ W 1,∞ on U i since F i is of class C 1 ∩ W 2,∞ on U i , det DF i 0 onŨ i , and hence the denominators in A i (x) are bounded away from zero on U i because of U i ⊂⊂Ũ i . Therefore, G i is of class C ∩ W 1,∞ on U i × B R for any R > 0.
The key idea is that, for any ( 0, since n(x) and ∇F i 3 (x) are parallel (and both non-zero). Thus, since the supports of the approximating functions ψ k shall be away from γ 0 T • and {0} × ∂Ω × R 3 , it is natural to consider the following C ∞ -function in the variables (t, G), that cuts off a region near the two sets where G 3 G 6 0 and where t G 3 0: For k ∈ N we then definẽ We should mention that, because of ζ i ∈ C ∞ c (U i ), i 0, . . . , m, the i-th summand is (by definition) zero if x U i . Note that we can apply the chain rule for η G i k since η k is smooth and G i ∈ W 1,1 (U i × B R ) for any R > 0. Therefore,ψ k is of class C ∩ W 1,∞ .
First we show that (4.3) holds forψ k (instead of ψ k ). By m i 0 ζ i 1 on Ω we have where C > 0 depends on the (finite) C 1 b -norms of ψ (and ζ i ) and where R > 0 is chosen such that ψ vanishes for t ≥ R or |v| ≥ R. For fixed i ∈ {1, . . . , m} and (t, x, v) ∈ R × U i × R 3 there hold the implications Therefore we have, recalling that 0 ≤ η ≤ 1, for l large. But then η G i k (t l , x l , v l ) 0 and therefore by (4.6) which is a contradiction. As for {0} × ∂Ω × R 3 , the proof works completely analogously.
There only remains one problem: The approximating functions are only of class C ∩ W 1,∞ with compact support and not of class C ∞ as desired (which corresponds to the fact that ∂Ω is only of class C 1 ∩ W 2,∞ and not necessarily smooth). To this end, take a Friedrich's mollifier ω ∈ C ∞ c R 7 , supp ω ⊂ B 1 , ∫ R 7 ω dvdxdt 1, and denote ω δ : δ −7 ω · δ for δ > 0.
Byψ k ∈ H 1 R 7 , we know that ω δ * ψ k converges toψ k for δ → 0 in H 1 R 7 . Moreover, since properties also hold for ω δ * ψ k instead ofψ k if δ is small enough. Choose 0 < δ k ≤ 1 such small and such that By p < 2, this implies where C > 0 depends on p, Ω, and R. After combining this with (4.5), noting thatψ k and ψ vanish for t ≥ R or |v| ≥ R and ω δ k * ψ k for t ≥ R + 1 (which implies the existence of r as asserted) or |v| ≥ R + 1, and setting we are finally done.
With this lemma, we can extend (2.1) to test functions ψ whose supports do not necessarily have to be away from γ 0 T • and {0} × ∂Ω × R 3 under a condition on the integrability of the solution.

Lemma 4.2. For fixed
Proof. Let 1 ≤ p < 2 satisfy 1 p + 1 q 1. In accordance with Lemma 4.1, let ψ k ⊂ Ψ T • approximate ψ with respect to the W 1,p t 2 x 1 v -norm, 0 < r < ∞ such that ψ and all ψ k vanish for t ≥ r, and define R : min{r, T • }. By assumption, (2.1) holds for ψ k for all k ∈ N. Hence there remains to show that we can pass to the limit k → ∞ in (2.1). First, we have for k → ∞. Note that this was the crucial estimate, for which we essentially needed the convergence of ψ k to ψ in the W 1,p t 2 x 1 v -norm. As for the boundary terms on γ ± T • , we first have for any t ∈ I T • , v ∈ R 3 , since Ω is bounded and ∂Ω of class C 1 . Therefore by n( for k → ∞, and the proof is complete. The next step is to show that (2.1) still holds if ψ does not depend on v. This is done via a cut-off procedure in v. Note that in the following lemma it is essential that f α is of class L 1 ∩ L 2 αkin locally in time.
(ii) If, additionally to the given assumptions, f α Proof. The proof works similarly to the proof of [11,Lemma 4.2.]. First, consider a test function ψ that may have support on ∂Ω.
Therefore, (2.1) holds for ψ m by Lemma 4.2. Now we can show that we may pass to the limit m → ∞ in all terms of (2.1) but the terms including integrals over γ ± T • . Let R > 0 such that ψ vanishes for t ≥ R. First, for (t, x, v) ∈ I T • × Ω × R 3 , we get the following estimate, which is again the crucial one: for m → ∞, since the last integral converges to zero by f α ∈ L 2 αkin [0, R] × Ω × R 3 . As for the term including the initial data, we see that The redundant divergence equations and the charge balance for m → ∞ as well by dominated convergence andf α ∈ L 1 Ω × R 3 . Now if supp ψ ⊂ [0, T • [ × R 3 \ ∂Ω , then ψ m vanishes on ∂Ω, too, and there vanish the integrals over γ ± T • for ψ m appearing in (2.1). Hence, (4.7) is satisfied. If the additional assumptions of (ii) hold, but ψ may not vanish on ∂Ω, we consider the integrals over γ ± T • : by dominated convergence and f α Therefore we obtain (4.8).
In the following, we denote

The redundant divergence equations and the charge balance
Proof. As for (i) and (ii), simply multiply (4.7) and (4.8) with e α and sum over α. As for (iii), Then ψ ∈ C ∞ I T • × R 3 with supp ψ ⊂ [0, T • [× R 3 compact. Therefore, Lemma 4.3 (ii) yields, after summing over α, from which the assertion follows immediately.
We can finally show the remaining parts of Theorem 2.5 with the help of Lemma 4.3; the redundancy of div x µH 0 has already been proved. To this end, assume Condition 2.4.
To prove (iii), let the additional assumptions stated there hold. The test function ϕ ∈ C ∞ c ]0, T • [ × R 3 may now not vanish on ∂Ω. Then we have ψ ∈ C ∞ I T • × R 3 with supp ψ ⊂ [0, T • [ × R 3 compact and Lemma 4.3 (ii) gives us, after multiplying with e α and summing over α, (4.13) We rewrite T ∂Ω ψ: Similarly as before, multiplying (4.11) and (4.13) with 4π and adding them to (4.10) yields Hence, div x (εE) 4π ρ int + ρ u + S ∂Ω on ]0, T • [ × R 3 in the sense of distributions. Hence, the f α and E, H of Theorem 2.3 satisfy these assumptions, and Theorem 2.5 (i), (ii) can be applied. However, the boundary values f α + constructed there only satisfy f α + ∈ L 1 lt γ + T • , dγ α for α 1, . . . , N ′ , i.e., the particles are subject to partially absorbing boundary conditions, and not necessarily for α N ′ + 1, . . . , N, i.e., the particles are subject to purely reflecting boundary conditions. Therefore, whether the statement of Theorem 2.5 (iii) is true for the solution of Theorem 2.3, remains as an open problem, unless N ′ N, i.e., all particles are subject to partially absorbing boundary conditions.
• The distribution S ∂Ω can be interpreted as follows: The terms j out ∂Ω (t, x) : where (t, x) ∈ I T • × ∂Ω, can be viewed as the outgoing and incoming boundary current density. Hence S ∂Ω can be rewritten as Thus, S ∂Ω measures how many particles have left and entered Ω up to time t. On the other hand, the distribution T ∂Ω measures how many particles leave and enter Ω at time t via We easily see that ∂ t S ∂Ω T ∂Ω on ]0, T • [ × R 3 in the sense of distributions, which corresponds to the fact that T ∂Ω appears as "a part of ∂ t ρ" in (4.9) and S ∂Ω appears as "a part of ρ" in (2.9).
• The global charge balance, see Corollary 4.4 (iii), can similarly been written as follows: for almost all t ∈ I T • .
• As mentioned in the introduction, in a more realistic model ε and µ should depend on f α , E, and H (maybe even nonlocally) and hence implicitly on time. In this situation, the weak formulation is the same as before, which is stated in Definition 2.1. If we assume ε, µ ∈ L ∞ loc I T • × R 3 ; R 3×3 (and suitably introduce initial values for ε, µ), viewed as explicit functions of t and x, the proofs of Theorem 2.5 and the lemmas before are still valid, and Theorem 2.5 remains true.
• Lastly, we emphasize that the results of this section hold, under the respective assumptions, for all weak solutions of (VM) in the sense of Definition 2.1 and not only for the solutions of Theorem 2.3.