Vlasov-Maxwell system with physical reasonableness constraints

Complete mathematical description on a physics system should contain not only equations reflecting physics laws but also some physical reasonableness constraints. For example, plasma is acknowledged to be described by Vlasov-Maxwell (V-M) system which contains 5 partial differential equations (PDEs): a Vlasov equation (VE) of the probability distribution function (PDF) f and 4 Maxwell equations (MEs). According to acknowledged rule, there should be $f\geqslant 0$ always. Therefore, strictly speaking, a complete mathematical description on plasmas contains 5 PDEs and an inequality (IE) f ≥ 0

$$\begin{eqnarray}\begin{array}{rcl}0 & = & \widehat{L}f\\ & = & \left[{\partial }_{t}+\upsilon \cdot {\rm{\nabla }}-\left[e\left[E\left(r,t\right)+\upsilon \times B\left(r,t\right)\right]\right]\cdot {\partial }_{p}\right]f\\ & = & \left[{\partial }_{t}+u\cdot {\rm{\nabla }}-e\left[E+u\times B\right]\cdot {\partial }_{p}\right]f\\ & & +\left[\upsilon -u\right]\cdot \left[{\rm{\nabla }}-{eB}\times {\partial }_{p}\right]f;\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\partial }_{t}E=e\int \upsilon {{fd}}^{3}\upsilon +{\rm{\nabla }}\times B;\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot E=-e\int {{fd}}^{3}\upsilon +{{ZeN}}_{i};\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot B=0;\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\times E=-{\partial }_{t}B,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&f\geqslant 0\forall r,p,t.\end{eqnarray} \tag{ 6 }$$

where p is the mono-variable function of υ: $p\left({var}\right)=\tfrac{{var}}{\sqrt{1-{var}\cdot {var}}}$, and N_i and Z are ionic density and charge respectively. In this standard model, the VE, equation (1), is formally a 1st-order partial differential equation (PDE) defined over a 7-D space $\left(t,x,y,z,{\upsilon }_{x},{\upsilon }_{y},{\upsilon }_{z}\right),4$ MEs are formally 1st-order PDEs defined over a 4-D space $\left(t,x,y,z\right)$, the IE is also defined over the 7-D space.

It is not trivial to emphasize these constraints. The IE (6) is not automatically valid. We introduce a mono-variable function ${f}_{\min }\left(t\right)\equiv \min \{f\left(r,p,\tau \right){| }_{\tau =t}\}$ to denote the minimum-value of f at each time point. According to the VE equation (1), for $f{| }_{t=0}\in {C}^{1}\left({R}^{3}\times {R}^{3},\left[{f}_{\min }{| }_{t=0},\infty \right)\right)$, there are $0={\rm{\nabla }}f{| }_{f={f}_{\min },t=0}={\partial }_{\upsilon }f{| }_{f={f}_{\min },t=0}$, which leads to, according to the VE equation (1), $0={\partial }_{t}f{| }_{f={f}_{\min }}$. Namely, for $f{| }_{t=0}\in {C}^{1}\left({R}^{3}\times {R}^{3},\left[{f}_{\min }{| }_{t=0},\infty \right)\right)$, if there is ${f}_{\min }{| }_{t=0}\gt 0$, there will be f > 0∀ r, p, t. But in many situations, $f{| }_{t=0}\notin {C}^{1}\left({R}^{3}\times {R}^{3},\left[{f}_{\min }{| }_{t=0},\infty \right)\right)$. For example, $f{| }_{t=0}$ is often chosen by researchers to be of a ramped density profile linear dependent on the space distance $| r-0|$: $\int f{| }_{t=0}{d}^{3}\upsilon =\beta | r|$. Thus, there will be

$$\begin{eqnarray*}&&{f}_{\min }{| }_{t=0}=0;\Longrightarrow {\rm{\nabla }}f{| }_{f={f}_{\min },t=0}\ne 0;\Longrightarrow {\partial }_{t}f{| }_{f={f}_{\min },t=0}\lt 0\exists \upsilon ;\Longrightarrow f{| }_{f={f}_{\min },t\gt 0}\lt 0\exists \upsilon .\end{eqnarray*}$$

Moreover, if f meets the VE, any mono-variable, but not limited to be non-negative-valued, function whose variable is f, also meets the VE. For example, sin(f). Therefore, not all initial condition $f\left(r,t,0\right)$ can warrant subsequently evolving f to meet the IE (6).

These examples well illustrate the necessity of basing investigation on plasmas on a complete mathematical description equations (1)–(6). Theoretical investigations on the V-M system (without the constraint) can be divided into two classes. One class is to apply analytic theories of PDE to the VE equation (1), which is taken as a 1st-order PDE defined over a 7-D space $\left(R\times {R}^{3}\times {R}^{3}\right)$, and study uniqueness, existence, regularity and conservation of various solutions with specialized global mathematical properties [1–7]. The other class is to apply standard computational mathematical approaches for numerically solving equations (1)–(5) [8–21]. The later class is more popular because it can give detailed expression of solutions in terms of space-time coordinates, which is necessary for theoretical-experimental comparison.

By now, all theoretical investigations are extension of differential equation theory, i.e. applying universal differential equation theoretical methods to PDEs reflecting physical laws. But the constraint IE reflecting physical reasonableness requirement is nearly untouched. Such an investigation style is logically sound: because PDEs with the constraint is a subset of PDEs, universal results valid for PDEs is also valid for PDEs with the constraint. But such universal results are too abstract to supply detailed solutions interested by physics community.

Such an investigation style is in a forward engineering manner: i.e., after obtaining universal results valid for PDEs and then picking out special results valid for PDEs with the constraint. Clearly, by now, such an forward engineering manner is still on the way (because the picking-out part is nearly not studied).

When there are multiple routes to achieve a same goal, the route with the highest efficiency is welcome. For studying such a complete mathematical description, beside above-mentioned forward engineering manner, there is a reverse engineering manner. That is, from the constraint IE, one can find allowed contour shapes which can warrant the IE to be fulfilled, and then based on these allowed contour shapes to study these PDEs. Compared with the forward engineering manner, the reverse engineering manner efficiently takes into account the constraint IE and targeted guides people attention being focused to some special results of PDEs. This makes finding detailed solutions interested by physics community to be in a short-cut, rather than in a long way.

Beside the VE, other equations of probability density function in many branches of applied sciences [22–25] should be also supplemented by a similar physical reasonable constraint inequality. Many authors [22–25] have noticed the importance of the non-negative-value constraint on the solved probability density functions and found a lot of solutions satisfying the constraint. Due to the difference between detailed forms of equations, those non-negative solutions cannot be transplanted to the V-M system, and hence efficiently finding non-negative solutions of f from the complete mathematical description on plasmas is worthy to be explored.

The purpose of this work is to display a strict, universal method, which is in such a reverse engineering manner, on equations (1)–(6). It is based on power series solution, and treats the VE equation (1) as a 'recurrence formula' relating many functions defined over $\left(r,t\right)$-space. These functions are expansion coefficient functions in the power series expression

$$\begin{eqnarray}&&f=\displaystyle \sum _{i\geqslant 0}{b}_{i}{\left[p\left(\upsilon \right)-p\left(\displaystyle \frac{\int \upsilon {{fd}}^{3}\upsilon }{\int {{fd}}^{3}\upsilon }\right)\right]}^{i}.\end{eqnarray} \tag{ 7 }$$

Especially, this method can, in the concrete, take the IE (6) into account.

An appropriate numerical experiment scheme on the complete mathematical description needs to re-consider the role of the VE and a more efficient expression of f. Strict theory has proven that any strict solution of the VE should have a spread over the υ-space [26]. This implies that f should be described by a set of functions $\{{b}_{i};0\leqslant i\leqslant \infty \}$ through power series [27–29]

$$\begin{eqnarray}&&f={b}_{0}\ast H\left(-| \upsilon -u| \right)+\displaystyle \sum _{i\gt 0}{b}_{i}{\{\left[p\left(\upsilon \right)-p\left(u\right)\right]\ast \exp \left(-{\left[p\left(\upsilon \right)-p\left(u\right)\right]}^{2}\right)\}}^{i},\end{eqnarray} \tag{ 8 }$$

where H is Heaviside function $H\left(x\geqslant 0\right)=1$ and $H\left(x\lt 0\right)=0$ and $\{{b}_{i\geqslant 0}\}$ are continuous functions of $\left(r,t\right)$. Thus, f is described by many 4-D functions, including $\int {{fd}}^{3}\upsilon$ and $\tfrac{\int \upsilon {{fd}}^{3}\upsilon }{\int {{fd}}^{3}\upsilon }$, those $\{{b}_{i\geqslant 1}\}$ reflect how the f deviates from a Dirac-type distribution with respect to υ. Thus, one can view the VE as playing a role of recurrence formula relating these expansion coefficient functions $\{{b}_{i\geqslant 0}\}$. Namely, this 1st-order 7-D PDE is equivalent to a set of 1st-order 4-D PDEs.

If further taking into account the mandatary constraint IE (6), one can find [28]

$$\begin{eqnarray}&&{b}_{1}\left(r,t\right)\equiv 0\end{eqnarray} \tag{ 9 }$$

to be mandatary (otherwise, there will be f < 0 regime nearby the super surface υ − u = 0 in the 6-D phase space). Thus, substituting the power series expression equation (8) into the VE and comparing terms proportional to ${\left[p\left(\upsilon \right)-p\left(u\right)\right]}^{i}$ order-by-order, one can find that for terms proportional to ${\left[p\left(\upsilon \right)-p\left(u\right)\right]}^{0}$, there should be

$$\begin{eqnarray}&&{\partial }_{t}{b}_{0}+u\cdot {\rm{\nabla }}{b}_{0}=0.\end{eqnarray} \tag{ 10 }$$

According to strict theory [26–29], equation (9) implies

$$\begin{eqnarray}&&\left[{\partial }_{t}p\left(u\right)+{eE}+{eu}\times B\right]=0;\end{eqnarray} \tag{ 11 }$$

Thus, the V-M system exhibits a universal architecture: it is expressed by 1st-order PDEs in infinite number. In the PDEs corresponded by the VE, one couples with MEs

$$\begin{eqnarray}&&{\partial }_{t}E={enu}+{\rm{\nabla }}\times B;\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot E=-{en}+{{ZeN}}_{i};\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot B=0;\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\times E=-{\partial }_{t}B,\end{eqnarray} \tag{ 15 }$$

to form a closed equation set of $\left(E,B,u\right)$, other PDEs in infinite number are expressions of those {b_i ≥ 0} in terms of $\left(E,B,u\right)$. Namely, one can determine 3 functions from equations (11)–(15), and then determine those {b_i ≥ 0} from above-mentioned other PDEs. Then, constructing a f through these solved functions according to equation (8), one can easily prove such a constructed f satisfying the VE. These other PDEs imply each b_i ≥ 2 to be a function of $\left(E,B,u,{b}_{0}\right)$, and ${b}_{0}{| }_{t=0}$ is determined through

$$\begin{eqnarray}&&{b}_{0}{| }_{t=0}={\left\{\int f\ast \delta \left(\upsilon -u\right){d}^{3}\upsilon \right\}}_{t=0}.\end{eqnarray} \tag{ 16 }$$

Namely, for f waiting to be solved, one can determine its two main features $\int {{fd}}^{3}\upsilon$, or ∇ · E, and u through equations (11)–(15) then determine other features, such as $\int f\ast \delta \left(\upsilon -u\right){d}^{3}\upsilon$, from the VE. Thus, roles of the VE and the MEs are effectively separated to avoid conflict. Clearly, the architecture displays that few fluid functions govern a few fluid functions. This agrees with basic mathematical principle. Opposite situation where a few govern few corresponds to the number of equations being more than the unknown and hence no solution exists.

The constraint IE (6) has two direct deductions: one is equation (9) and the other is

$$\begin{eqnarray}&&\int {{fd}}^{3}\upsilon =n={{ZN}}_{i}-\displaystyle \frac{{\rm{\nabla }}\cdot E}{e}\geqslant 0\forall r,t.\end{eqnarray} \tag{ 17 }$$

Equations (11)–(15) is a 1-order PDE set of three functions $\left(E,B,u\right)$. It can also be expressed as the combination of the IE (17) and a 2-order PDE set of two functions $\left(B,u\right)$ by expressing E in terms of p and B through equation (11) and substituting it into MEs equations (12)–(15) Note that the constraint equation (17) becomes to be in terms of p and B

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\partial }_{t}p+{\rm{\nabla }}\cdot \left[u\times B\right]+{{ZN}}_{i}\geqslant 0\forall r,t.\end{eqnarray} \tag{ 18 }$$

Now, what we dealt with is MEs equations (12)–(15)in terms of two functions p and B, which is an 2-order PDE set of $\left(B,u\right)$, and the constraint IE (18).

The role of the constraint IE of determining the shape of space-time contours of related physical fields has been not fully appreciated and utilized yet. For such a partial differential IE (PDIE) like equation (18), we can use independent variables, similar to those in [22, 23], constructed from $\left(r,t\right)$

$$\begin{eqnarray}&&{\xi }_{+}\equiv \int h\left(r\right){dr}+\nu t;{\xi }_{-}\equiv \int h\left(r\right){dr}-\nu t,\end{eqnarray} \tag{ 19 }$$

where the function h(r) represents the space shape of the static ionic density profile N_i ≡ N_i,0h(r) (N_i,0 is a constant), and ν is a constant, to rewrite the PDIE as a general form (where c_{i, j} are constants)

$$\begin{eqnarray}&&\displaystyle \sum _{i,j}{c}_{i,j}{\partial }_{{\xi }_{+}^{i}}^{i}{\partial }_{{\xi }_{-}^{j}}^{j}\left[\displaystyle \frac{1}{2}\left({\partial }_{{\xi }_{+}}-{\partial }_{{\xi }_{-}}\right)p+u\times B\right]+{{ZN}}_{i,0}\geqslant 0\forall r,t.\end{eqnarray} \tag{ 20 }$$

If expressing p and u × B as polynomials of ξ₊ and ξ₋, we can easily find that the IE (18) has confined p and u × B to be mono-variable functions of a same independent variable ξ₊ or ξ₋, i.e.

$$\begin{eqnarray}&&p\equiv p\left({\xi }_{+}\right);B\equiv B\left({\xi }_{+}\right);\end{eqnarray} \tag{ 21 }$$

or

$$\begin{eqnarray}&&p\equiv p\left({\xi }_{-}\right);B\equiv B\left({\xi }_{-}\right).\end{eqnarray} \tag{ 22 }$$

Otherwise, the IE cannot be valid ∀r, t because the equality ${\sum }_{i,j}{c}_{i,j}{\xi }_{+}^{i}{\xi }_{-}^{j}+{{ZN}}_{i,0}=0$ will define a super curve on the ${\xi }_{+}---{\xi }_{-}$ plane and hence divide the plane into two parts: one for > 0 and the other for < 0. In contrast, if the lefthand side of the IE (20) is mono-variable function and thus the PDIE is indeed an ordinary differential IE (ODIE), one can easily find that if the minimum of this mono-variable function is above the constant $-{{ZN}}_{i,0}$, the IE is valid ∀r, t. Namely, the IE has determined the shape of the space-time contours of the p-field and B-field. Moreover, why traveling wave solutions of other equations of probability density function are emphasized [22–25] might be reasonably attributed to the IE' requirement.

For $h\left(r\right)\equiv 1$-case, equations (12)–(15) implies a 2nd-order PDE of W + u × B [27–30]

$$\begin{eqnarray}&&\{{\partial }_{{tt}}-{\rm{\Delta }}+{\partial }_{t}\left[u{\rm{\nabla }}\cdot \right]+{\rm{\nabla }}{\rm{\nabla }}\cdot \}\left[W+u\times B\right]=0,\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&W\equiv {\partial }_{t}p\left(u\right)+{{ZN}}_{i}r/3.\end{eqnarray} \tag{ 24 }$$

If denoting the angle between u and B as θ, we can find that W is θ-independent while u × B θ-dependent [30]. Because the allowed values of θ are over a broad regime, we thus obtain two same 2nd-order PDEs [27–30]

$$\begin{eqnarray}&&\{{\partial }_{{tt}}-{\rm{\Delta }}+{\partial }_{t}\left[u{\rm{\nabla }}\cdot \right]+{\rm{\nabla }}{\rm{\nabla }}\cdot \}W=0;\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\{{\partial }_{{tt}}-{\rm{\Delta }}+{\partial }_{t}\left[u{\rm{\nabla }}\cdot \right]+{\rm{\nabla }}{\rm{\nabla }}\cdot \}\left[u\times B\right]=0.\end{eqnarray} \tag{ 26 }$$

Applying ${\rm{\nabla }}\cdot$ to this equation set, we can find that ∇ · W and ${\rm{\nabla }}\cdot \left[u\times B\right]$ also meet a same 2nd-order PDE and hence this fact implies a linear relationship between them [29, 30]

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left[u\times B\right]={c}_{0}+{c}_{1}\ast {\rm{\nabla }}\cdot {\partial }_{t}p;\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\left[u\times B\right]={c}_{0}\ast \displaystyle \frac{r}{3}+{c}_{1}\ast {\partial }_{t}p,\end{eqnarray} \tag{ 28 }$$

where c₀ and c₁ are constants. Because the linear relation should hold $\forall t$ (i.e., including t = 0), the values of c₀ and c₁ can be determined by initial condition $\left({\partial }_{t}p,u,B\right){| }_{t=0}$ [30, 31]:

$$\begin{eqnarray}&&{c}_{1}=-(1+\alpha )\ast \max \{{\{{\rm{\nabla }}\cdot {\partial }_{t}p+{{ZN}}_{i}\}}_{t=0}\};\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{c}_{0}={{ZN}}_{i}{c}_{1},\end{eqnarray} \tag{ 30 }$$

where α is a positive constant and the function max represents the maximum of allowed values of its variable. Thus, the constraint IE (18) becomes

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\partial }_{t}p+{{ZN}}_{i}\leqslant 0\forall r,t.\end{eqnarray} \tag{ 31 }$$

As above-mentioned, the IE (18) determines these PDEs to be finally converted to ODEs with respect to the independent variable ξ₊ or ξ₋. Note that equation (25) is a 2nd-order PDE of W or a 3rd-order PDE of p. A time-conservation law can be derived from this 3rd-order PDE [30]. Thus, equations (12)–(15), (18) is converted into a 2nd-order PDE of u and expressions of E and B in term of u [28, 30] and the constraint IE (18) (for simplicity, we do not write out their corresponding ODEs here):

$$\begin{eqnarray}&&{\partial }_{{tt}}p+\left[\left({\partial }_{t}{\rm{\nabla }}\cdot p\right)\displaystyle \frac{p}{\sqrt{1+{p}^{2}}}\right]+{\rm{\nabla }}\times {\rm{\nabla }}\times p=-\left[{{uZN}}_{i}\right]+{POT};\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\partial }_{t}p+{{ZN}}_{i}r/3=\displaystyle \frac{1}{{c}_{1}}\left[u\times B\right];\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&-E={\partial }_{t}p+u\times B;\end{eqnarray} \tag{ 34 }$$

where POT is a vector determined by initial conditions of u, which is indeed determined by the initial condition of $\left(E,B\right)$. This is because equation (25) is a 3rd-order PDE of p. Initial conditions for determining solutions of equation (25) include 3 terms: $\left(u,{\partial }_{t}u,{\partial }_{{tt}}u\right){| }_{t=0}$, and two of them $\left({\partial }_{t}u,{\partial }_{{tt}}u\right){| }_{t=0}$ are affected by $\left(E,B\right){| }_{t=0}$ through equation (11).

The positive-valued constant α is also determined by initial condition. The vector POT can be viewed as a constant if initially electromagnetic (EM) fields is of nearly uniform transverse intensity distribution. If initially the fields is a pure EM fields such as light, there will be $\alpha =\tfrac{{\omega }_{B}}{{\omega }_{p}}\tfrac{{\omega }_{B}}{\omega }\tfrac{1}{\sqrt{1+{\left(\tfrac{{\omega }_{B}}{\omega }\right)}^{2}}}-1$ and ${POT}=\max \left\{{\partial }_{{zz}}A{| }_{t=0}\right\}$ [30, 31], where A, ω, ω_B are the vector potential, the frequency and the cyclotron frequency of the initial EM fields and ω_p is the plasma frequency. For underdense plasmas, ω_p < ω can warrant $\alpha \gt \displaystyle \frac{{\omega }_{B}}{\omega }\displaystyle \frac{{\omega }_{B}}{\omega }\displaystyle \frac{1}{\sqrt{1+{\left(\tfrac{{\omega }_{B}}{\omega }\right)}^{2}}}-1=\sqrt{1+{\left(\tfrac{{\omega }_{B}}{\omega }\right)}^{2}}-\tfrac{1}{\sqrt{1+{\left(\tfrac{{\omega }_{B}}{\omega }\right)}^{2}}}\gt 0$.

Complete mathematics description on plasmas is the V-M system constrained by an IE f ≥ 0∀r, υ, t . This is a model containing a 1st-order PDE and an IE defined over a 7-D space, four 1st-order PDEs over a 4-D space. Through the power series solution method, the constraint IE f ≥ 0∀ r, υ, t can be embodied through a mandatary fluid equation that expansion coefficient functions should satisfy and a macroscopic constraint IE n ≥ 0∀ r, t. The 7-D PDE plays a role of 'recurrence formula' relating these expansion coefficient functions (defined over the 4-D space). As a consequence, the complete mathematics description displays a universal architecture which refers to 3 fields functions defined over the 4-D space governing these expansion coefficient functions in infinite number. What 3 governing fields satisfy also contains the macroscopic constraint IE n ≥ 0∀ r, t. A strict method of solving fluid equations of 3 governing fields and the IE n ≥ 0∀ r, t is found. The IE n ≥ 0∀ r, t determines/affects the shape, or geometric characteristics, of contours of the space-time profile fluid velocity field.

In short, with increasing applications of statistics theory/viewpoint to various science disciplines, beside the detailed form of the evolution equation of the probability density function, (which might vary with respect to studied topics), non-negative-value constraint on the probability density function is mandatory to be satisfied. Beside directly seeking for non-negative solution of the probability density function [22–25], more in-depth treatment on the constraint is desired. Compared with the situation [22] where the probability density function is only space-time dependent, or a 4-D function, and its evolution equation is a PDE, the probability density function in plasma physics, a 7-D function, is more complicated and its evolution equation is indeed a differential-integral equation (because of MEs). Especially, the evolution equation involves in the definite integral over an independent variable (such as υ). In such a class of situations, directly seeking for non-negative solutions of the evolution equation is not an easy task, and the dependence of the probability density function on the integrated independent variables needs in-depth consideration. Therefore, it is necessary to analyze what the constraint implies and derive more concrete formulas from the constraint. This work displays an example of how to derive more concrete formulas from the constraint. Obviously, these concrete formulas are favorable to efficiently find non-negative solutions.