The existence of numerous, up to astronomical figures, interacting charged particles adds difficulty to the simulation on plasma. For Vlasov-Maxwell (V-M) description on plasma   , there are usually two schemes of simulation. One is well-known fluid scheme  in which the microscopic distribution function
is represented by a series of moments obeying a series of fluid equations , where is Vlasov equation (VE)    , the operator , represents the Lorentz force, and . At present, the fluid scheme is limited by the difficulty in handling such a series whose members are in infinite
number. People often resort to truncation approximation to the series and hence lower down the scientific validity of results derived from the fluid scheme. Another scheme is to express as a mono-variable function of the invariants of single particle motion  -  . Because the invariants of single particle motion reflect a relation among u and the self-consistent fields and , to know final expression of f in term of , we need exact solutions of E and B. According to Maxwell equations (MEs), E and B are coupled with . Even an exact expression of f in term of the invariants of single particle motion is available, it is less helpful for obtaining exact solutions of E and B because the integral of f over u-space will lead to a space-time function whose expression in term of E and B is very complicated. If still trying to solve exact solutions of E and B along this way, approximation on the integral of f over u-space is inevitable.
Above-mentioned difficulty in studying plasma on the basis of the V-M description causes people to try Newton-Maxwell (N-M) description in which each realistic particle meets a relativistic Newton equation (RNE) . Because realistic particles studied are usually in astronomical figures, for making simulation to be practical, people often merge numerous realistic particles into a macroparticle and hence the number of macroparticles studied is smaller several magnitudes than that of realistic particles   . Because the merging is according to initial positions, regardless of initial velocities, of realistic particles, this automatically implies a rigid-macroparticle approximation which means each macroparticle keeping its realistic members always as an entirety.
Clearly, the rigid-macroparticle approximation corresponds to a distorted picture. True picture is that each macroparticle, due to differences among velocities of its composite realistic particles, will break into pieces of different destinations at next time point, rather than moving as an entirety to a destination at next time point. Such a distorted picture sheds an uncertainty on the scientific validity of its yielding numerical results.
The purpose of this work is to solve such an inherent drawback of the N-M description on plasma. We will display in details that it is feasible to strictly simulate plasma beyond the rigid-macroparticle approximation. Such a feasibility arises from an inherent agreement between the V-M description and the N-M one, which is exhibited by a same equation reflecting the relation between
2. Materials and Methods
We start from the V-M description, which consists of Vlasov equation (VE) and MEs. The VE reflects that the space-time evolution of f is governed by the self-consistent fields . The MEs reveal the dependence of on f. According to strict theory, exact solutions of a VE always has a thermal spread or a spread over the u-space  . Because the VE can lead to a set of fluid
equations or from to , and in each i case, the term will involve in moments from to , this implies that all moments form an open set (because the number of its members is infinite).
According to MEs, couples with and formally decouples with all . According to an open fluid equations set , each equation reflects a universal relation among , and at least . This seems that exact solutions of all higher-order moments are necessary for that of . On the other hand, because and are nonlinear functions of u and hence , a term in the equation , is dependent on moments from to , each equation will imply a relation among moments from to and thus all equations from to will mean a matrix describing the relation among all moments from to . Clearly, obtaining exact solutions of from such an open equation set is impossible, and, as discussed below, also unnecessary.
Another open set , where , can be
defined naturally through the M-set   . Clearly, and automatically exist. Each equation can be expressed through the D-set
and coefficients are known functionals of . For example, and can lead to to be expressed through Equation (1). Starting from the case, we can formally obtain an expression of in all terms , and then substituting it into the case and formally obtain an expression of in all terms ,.... Finally, we will find that all are determined by and all coefficients . Namely, the open equation set does not lead to a substantial constraint on .
According to MEs, depend on and is independent of the D-set. Therefore, exact solutions of the D-set is not a necessary condition for those of . The open equation set is only responsible for relations among those and cannot has an effect on .
Strict mathematical theory have revealed that exact solutions of can be obtained through a functional of f,
and (see the appendix of Ref.  ). For simplicity of symbols, we denote as W and F as , where
. From the definition of F, there is always . In the appendix of Ref.  or    , strict mathematics
have proven . For convenience of readers, we make a clearer presentation of detailed proof as below:
It is easy to directly verify following relations
where we have used the relation , and
When deriving the Equation (3), we have used two relations:
and . Moreover, because
where , is a power series of and does not have zero-order term, it is easy to verify that there exists
This explains another relations used in deriving the Equation (3)
Thus, by shifting the term , which is equal to , from righthand side to lefthand side of the “=”, we can strictly rewrite  as
Clearly, it has a strict solution , or . Here, why we only choose the exact solution and ignore other exact solutions such as is explained as follows: if we express f as a power series of , i.e. , because u is
defined as , or , and must be satisfied,
there usually should be . Thus, if substituting this power series into the VE and comparing terms order-by-order, we will find that for terms , because of , there will be . According to this power series, is just . Consequently, because implies , this forces us to choose the exact solution .
Detailed expression of is
according to standard procedure, it can cause to be equivalent/reduced to
Equation (9) indicates that F cannot leads to a CE (because it does not satisfy the VE). Moreover, although those higher-order moments still formally meet a set of equations in infinite-number the Dirac function dependence of F on u make these equations in infinite-number indeed to be equivalent to a same equation, Equation (10). This can be easily verified by simple algebra.
Equation (10), , and the CE enable each equation in the open set
, to yield a relation among all higher-order off-center
moments . For example, can reflect a relation, in a general form Equation (1), among higher-order moments in infinite-number. Clearly, for solving those moments in infinite-number, merely an equation or Equation (1) at -case is insufficient and more similar equations in infinite-number, or Equation (1) at all -cases, are required. The open set can in principle yield the expression of every higher-order off-center moments in term of first two moments and , or in term of , which can be solved from Equation (10) and 4 MEs
Therefore, there is a theorem: For any V-M system, its meets . Namely, means . Namely, the f of any V-M system can be solved through two equations
where is a number. Due to and the Dirac function dependence of F on u, all will be equivalent to a same equation 
Due to the universal validity of Equations (12, 13) for any a within and the equality , Equations (12) & (13) will lead to a balance
Namely, the convective term combines with to determine those through Equation (1) at all i-cases, or an matrix equation of the raw vector . As previously pointed out, the equation reflect a relation among all . To determine these , the whole set is required. If simply discarding those -related terms, we will make the remained convective term to affect u and hence . This leads to inexact solutions of .
Equation (10) can also be derived from Newton-Maxwell (N-M) system. As proven in the appendix of Ref.  , for a group of electrons , we can always define two fields (in Lagrangian expression)
which are fluid velocity field and relative velocity field respectively. Especially, there is always a formula
where º means the formula is automatically valid for any . If applying to this automatically valid formula, we will obtain another automatically valid formula
No matter what the relative field is, the relativistic Newton equation (RNE) of every electron is always valid. Thus, the condition for the RNE being valid under arbitrary value of the RV-field, of course including a common value: , will lead to the Lagrangian expression of Equation (10).
This implies a solution to the dependence of calculated in macroparticle dynamics simulation   on the graininess parameter G, which reflects how many realistic particles are contained in a macroparticle. Namely, G-independent can be obtained through Equation (10) and MEs. Once this dependence is removed, smaller G-parameter will correspond to a reliable, finer description on the projection of f on u-space.
Many typical exact solutions of are well known. For example, an exact solution of Equations (10) & (11) with a constraint and a condition could describe light (or pure transverse electromagnetic wave) in vacuum: , where meets . Likewise, Equations (9) & (11) with a condition can describe light-matter interaction. Expressing E in term of u and B   , we find that Equations (10) & (11) will yield
where we have used the relation . As analyzed elsewhere    , Equation (20) can be further written as
Moreover, any charged particles system should have a constraint or because the density cannot be negative-valued. Equation (21) directly implies
where POT is a constant vector determined by the initial condition of the interaction.
There are some examples of the application of Equation (10) in calculating plasma self-consistent fields   . Some typical analytical formula of exact solutions of f can be found elsewhere  .
We have outlined, with strict mathematical proof, a feasible scheme of simulating plasma beyond rigid-macroparticle approximation. It enables exact solutions of the self-consistent fields to be available. Consequently, exact solutions of microscopic distribution functions f are warranted. The scheme is of a universal application value to plasma and beam physics.
Conflict of Interest
There is no conflict of Interest in this work.
This work is supported by the Natural Scientific Fund no 11374318.