Asymptotic-preserving neural networks for multiscale
  Vlasov-Poisson-Fokker-Planck system in the high-field regime
The Vlasov-Poisson-Fokker-Planck (VPFP) system is a fundamental model in
plasma physics that describes the Brownian motion of a large ensemble of
particles within a surrounding bath. Under the high-field scaling, both
collision and field are dominant. This paper introduces two
Asymptotic-Preserving Neural Network (APNN) methods within a physics-informed
neural network (PINN) framework for solving the VPFP system in the high-field
regime. These methods aim to overcome the computational challenges posed by
high dimensionality and multiple scales of the system. The first APNN method
leverages the micro-macro decomposition model of the original VPFP system,
while the second is based on the mass conservation law. Both methods ensure
that the loss function of the neural networks transitions naturally from the
kinetic model to the high-field limit model, thereby preserving the correct
asymptotic behavior. Through extensive numerical experiments, these APNN
methods demonstrate their effectiveness in solving multiscale and high
dimensional uncertain problems, as well as their broader applicability for
problems with long time duration and non-equilibrium initial data.
Shi Jin, Zheng Ma, Tian-ai Zhang
Asymptotic-preserving neural networks for multiscale Vlasov-Poisson-Fokker-Planck system in the high-field regime

###### Abstract

The Vlasov-Poisson-Fokker-Planck (VPFP) system is a fundamental model in plasma physics that describes the Brownian motion of a large ensemble of particles within a surrounding bath. Under the high-field scaling, both collision and field are dominant. This paper introduces two Asymptotic-Preserving Neural Network (APNN) methods within a physics-informed neural network (PINN) framework for solving the VPFP system in the high-field regime. These methods aim to overcome the computational challenges posed by high dimensionality and multiple scales of the system. The first APNN method leverages the micro-macro decomposition model of the original VPFP system, while the second is based on the mass conservation law. Both methods ensure that the loss function of the neural networks transitions naturally from the kinetic model to the high-field limit model, thereby preserving the correct asymptotic behavior. Through extensive numerical experiments, these APNN methods demonstrate their effectiveness in solving multiscale and high dimensional uncertain problems, as well as their broader applicability for problems with long time duration and non-equilibrium initial data.

## 1 Introduction

The Vlasov-Poisson-Fokker-Planck (VPFP) system is a fundamental kinetic model in plasma physics, which takes into account the interactions between electrons and a surrounding bath via the Coulomb force [6]. It consists of a Liouville equation with a Fokker-Planck operator in velocity space, along with a Poisson equation corresponding to the electrostatic force [43]. The high-field limit of the kinetic equations characterizes a usual scenario where the external field, such as an electric field, is strong enough to balance the collision term [14]. Section 2 provides a detailed exposition of the VPFP system and its high-field limit.

The high dimensionality of the phase space in the kenitie equations (6 dimensions plus time for the general cases) raises a significant challenge in simulating the VPFP system. Moreover, the uncertainties of the kinetic models can introduce more dimensions [16, 37, 38, 18]. The uncertainties may arise from the collision kernels, scattering coefficients, initial or boundary data, source or forcing terms and so on, which are important for the application of the kinetic models to practical systems [2, 3, 15, 22]. Many efforts have been devoted to tackling the curse of dimensionality. Classical techniques, such as the Monte Carlo
methods [13, 31], are hindered by their low-order accuracy or slow convergence rates. Recently, machine learning, especially deep neural networks (DNNs), shows potential in resolving the Partial Differential Equations (PDEs) of high dimensionality [21, 23, 30, 24, 33]. The key idea is to approximate the objective functions by DNNs, and optimize the network parameters by minimizing the loss function encoding the PDEs and initial-boundary conditions. Many work has been devoted to building suitable loss functions to develope DNN methods. The popular methods are the Physics-Informed Neural Networks (PINNs) [40] and the Deep Galerkin Method (DGM) [41], which minimize the \(L^{2}\)-residual of the PDEs and the initial-boundary conditions. The Deep Ritz Method [47] was designed to solve PDEs with variational structures by leveraging its Ritz formulation. Moreover, the new area of operator learning has risen, which involves developing novel learning frameworks to approximate the solution operator of PDEs. Such techniques include deep operator networks (DeepONets) [7, 33], neural operators [23, 28] and associated variants tailored for physics problems [29, 44]. For supplementary information on machine learning methods for solving PDEs, interested readers are advised to consult the comprehensive review [45]. Applications of these methods to the kinetic equations are detailed in references [8, 11, 27, 29, 32].

Besides the curse of dimensionality, the multiscale nature of the VPFP system with high-field limit poses another challenge. General DNN methods, like PINNs, have difficulties in dealing with multiscale problems [17]. Specifically, PINNs only capture the leading order term related to a single scale in loss functions when the scale parameter is small [20], thus fail to capture the correct asymptotic limit equation of the VPFP system in the high-field regime as shown in Section 3. Asymptotic-Preserving (AP) schemes are effective numerical schemes for multiscale problems. AP schemes preserve the continuous asymptotic limit in a numerically uniformly stable way [15]. In recent years, AP schemes are widely adopted for kinetic and hyperbolic equations with multiple time and space scales [14], and also applied in the plasma [10] and the high-field regime [4, 9, 19, 48]. However, the classical AP schemes still encounter challenges in the high dimensional cases.

To address both the curse of dimensionality and multiscale difficulty, the primary goal of this work is to develope DNN methods with AP property, which are applied to VPFP system in the high-field regime. We propose two Asymptotic-Preserving Neural Network (APNN) methods for the VPFP system. The first is based on the micro-macro model of the VPFP system [9]. By reformulating the loss function of the vanilla PINN method based on the micro-macro decomposition, we construct an associated loss function with AP property. This method has been successfully applied to the linear transport equations [20], the steady radiative transfer equations [34] and so on [17, 26]. However, this method necessitates initial data close to the local Maxwellian, since its foundation is the micro-macro decomposition. Since kinetic models are more useful when the system is far from the equilibrium initial data [19], we try to retain this advantage. Hence, we propose the second APNN method, which is free from the explicit expression of the local Maxwellian and the requirement for the equilibrium initial data. This method is based on the mass conservation law, and has the potential for a broader applicability for the cases with long time duration or non-equilibrium initial data.

Although our numerical experiments are conducted in one-dimensional physical and velocity spaces, we do include the examples that involve high-dimensional uncertainties to demonstrate the ability of our APNNS for high-dimensional problems.

The paper is organized as follows: Section 2 introduces the VPFP system and its high-field limit, along with the micro-macro model of the VPFP system. Section 3 presents our two proposed APNN methods and demonstrates their AP property. In Section 4, we evaluate and compare the performance of the two proposed APNN methods through a series of numerical examples. We conclude the paper in Section 5.
The VPFP system and its high-field limit

This section introduces the Vlasov-Poisson-Fokker-Planck (VPFP) system and its high-field limit. The micro-macro model for the VPFP system is also presented, which provides a suitable framework for deriving the high-field limit.

### The VPFP system

In kinetic equations, a strong external field, such as the electric or gravitational field, often exists and counteracts the collision term, leading to the high-field limit [5]. In electrostatic plasmas for example, the electrons interact with the surrounding bath through the Coulomb force. Thus the time evolution of the electron distribution function \(f:(t,x,v)\in\mathbb{R}_{+}\times\mathbb{R}^{N}\times\mathbb{R}^{N}\to \mathbb{R}_{+}\) is governed by the VPFP equations, with the action of a self-consistent potential \(\phi(t,x)\). Specifically, the VPFP system is:

\[\partial_{t}f+v\cdot\nabla_{x}f-\frac{1}{\varepsilon}\nabla_{x}\phi\cdot\nabla _{v}f=\frac{1}{\varepsilon}\nabla_{v}\cdot[vf+\nabla_{v}f]:=\frac{1}{ \varepsilon}\mathcal{Q}(f), \tag{2.1a}\] \[-\triangle_{x}\phi(t,x)=\rho(t,x)-h(x), \tag{2.1b}\]

where

\[\rho(t,x)=\int_{\mathbb{R}^{N}}f(t,x,v)\mathrm{d}v \tag{2.2}\]

is the density of electrons. \(\mathcal{Q}\) is defined as linear Fokker-Planck operator:

\[\mathcal{Q}(f)(t,x,v)=\nabla_{v}\cdot\left[vf(t,x,v)+\nabla_{v}f(t,x,v)\right]. \tag{2.3}\]

The function \(h(x)\) represents a given positive background charge density. Hence the assumed global neutrality relation is:

\[\int_{\mathbb{R}^{N}}\int_{\mathbb{R}^{N}}f^{0}(x,v)\mathrm{d}x\;\mathrm{d}v =\int_{\mathbb{R}^{N}}h(x)\mathrm{d}x. \tag{2.4}\]

The parameter \(\varepsilon\) denotes the ratio between the mean free path and the Debye length. The limiting process \(\varepsilon\to 0\) indicates the high-field limit of the VPFP system. In the high-field regime, the strong forcing term balances the Fokker-Planck diffusion term [1], which is different from the low-field limit (or named as parabolic limit) [39].

### The high-field limit

In this section, we give the limit equation of the VPFP system 2.1a-2.1b at \(\varepsilon\to 0\). First, integrating the Vlasov equation 2.1a over \(v\) in \(\mathbb{R}^{N}\) gives:

\[\partial_{t}\int_{\mathbb{R}^{N}}fdv+\nabla_{x}\cdot\int_{\mathbb{R}^{N}}vfdv -\frac{1}{\varepsilon}\int_{\mathbb{R}^{N}}\nabla_{v}\cdot(\nabla_{x}\phi fdv )=\frac{1}{\varepsilon}\int_{\mathbb{R}^{N}}\nabla_{v}\cdot(vf+\nabla_{v}f) \,dv. \tag{2.5}\]

Integrating by parts, one has:

\[\partial_{t}\rho+\nabla_{x}\cdot j=0, \tag{2.6}\]

where the flux \(j\) is defined as:

\[j=\int_{\mathbb{R}^{N}}vf(t,x,v)\mathrm{d}v. \tag{2.7}\]
Next, one multiplies the Vlasov equation (2.1a) by \(v\) then integrates over \(\mathbb{R}^{N}\), and takes the limit \(\varepsilon\to 0\), yielding:

\[0=\int_{\mathbb{R}^{N}}f\nabla_{x}\phi+vf+\nabla_{v}fdv. \tag{2.8}\]

Deriving this equation gives:

\[j=-\rho\left(\nabla_{x}\phi\right). \tag{2.9}\]

Substituting it in equation (2.6) yields the high-field limit equation:

\[\left\{\begin{array}{l}\partial_{t}\rho-\nabla_{x}\cdot(\rho\nabla_{x}\phi)= 0,\\ -\Delta_{x}\phi=\rho-h(x).\end{array}\right. \tag{2.10}\]

The rigorous proof for the high-field limit of the VPFP system in one-dimension can be found in [12, 36].

### The micro-macro model and the high-field limit

This section presents the micro-macro model, from which the high-field limit can be readily derived. First, an equivalent formulation of the VPFP system (2.1a)-(2.1b) is:

\[\partial_{t}f+v\cdot\nabla_{x}f =\frac{1}{\varepsilon}\nabla_{v}\cdot[(v+\nabla_{x}\phi)f+ \nabla_{v}f]:=\frac{1}{\varepsilon}\mathscr{L}f, \tag{2.11}\] \[-\Delta_{x}\phi =\rho-h(x).\]

The linear operator \(\mathscr{L}\) retains the formulation of a Fokker-Planck operator dependent on \(\nabla_{x}\phi\). The null space of \(\mathscr{L}\) is:

\[\mathcal{N}(\mathscr{L})=\text{Span}\{\mathscr{M}\}=\{f=\rho\mathscr{M}, \text{where}\,\rho:=\langle f\rangle\}, \tag{2.12}\]

where \(\mathscr{M}\) is the shifted Maxwellian depending on \((t,x)\) through the potential \(\phi\) in the form:

\[\mathscr{M}(t,x,v)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{|v+\nabla_{x}\phi(t, x)|^{2}}{2}\right). \tag{2.13}\]

Notation \(\langle f\rangle\) is defined as:

\[\langle f\rangle=\int_{\mathbb{R}^{N}}f(v)dv. \tag{2.14}\]

The rank of \(\mathscr{L}\) is:

\[\mathscr{R}(\mathscr{L})=(\mathcal{N})^{\perp}(\mathscr{L})=\{f\text{ such that }\langle f\rangle=\int_{\mathbb{R}^{N}}fdv=0\}. \tag{2.15}\]

The idea for the micro-macro model begins with the decomposition of \(f\) into the equilibrium part and the non-equilibrium part [9, 25]:

\[f=\rho\mathscr{M}+\varepsilon g, \tag{2.16}\]

where \(\rho\mathscr{M}\) is the equilibrium part of \(f\). The non-equilibrium part \(g\) satisfies \(\langle g\rangle=0\), which implies \(g\in\mathscr{R}(\mathscr{L})\). The orthogonal projector \(\Pi\) in \(L^{2}\left(\mathscr{M}^{-1}dv\right)\) onto \(\mathcal{N}(\mathscr{L})\) is defined by:

\[\Pi\varphi=\langle\varphi\rangle\mathscr{M}. \tag{2.17}\]

In deriving the micro-macro model, we first consider the one-dimensional VPFP system formulated as:

\[\partial_{t}f+v\partial_{x}f=\frac{1}{\varepsilon}\partial_{v}\left[(v+ \partial_{x}\phi)f+\partial_{v}f\right]=:\frac{1}{\varepsilon}\mathscr{L}f. \tag{2.18}\]
Applying the micro-macro decomposition (2.16) to the equation (2.18) yields:

\[\partial_{t}(\rho\mathcal{M})+\varepsilon\partial_{t}g+v\partial_{x}(\rho \mathcal{M})+\varepsilon v\partial_{x}g=\mathcal{L}g. \tag{2.19}\]

Operating \((I-\Pi)\) on equation (2.19) gives:

\[(I-\Pi)\left[\partial_{t}(\rho\mathcal{M})+v\partial_{x}(\rho\mathcal{M}) \right]+\varepsilon(I-\Pi)\left[\partial_{t}g+v\partial_{x}g\right]=(I-\Pi) \mathcal{L}g.\]

Since \(\Pi(g)=\Pi\left(\partial_{t}g\right)=\Pi(\mathcal{L}g)=0\)[9], and

\[\Pi\left[\partial_{t}(\rho\mathcal{M})\right]=\Pi\left[\mathcal{M}\partial_{ t}\rho-\rho\left(\partial_{xt}\phi\right)(v+\partial_{x}\phi)\mathcal{M} \right]=\mathcal{M}\partial_{t}\rho,\]

one gets

\[\varepsilon\partial_{t}g+\varepsilon(I-\Pi)\left(v\partial_{x}g\right)= \mathcal{L}g-(I-\Pi)\left(v\partial_{x}(\rho\mathcal{M})\right)-\rho( \partial_{t}\mathcal{M}). \tag{2.20}\]

Then the \(\Pi\) projection (2.17) is applied to the equation (2.19)

\[\partial_{t}\rho+\partial_{x}(v(\rho\mathcal{M}))+\varepsilon\partial_{x}(vg) =0. \tag{2.21}\]

Therefore, one reformulates the original VPFP system (2.1a)-(2.1b) as:

\[\begin{cases}\partial_{t}g+(v\partial_{x}g-\partial_{x}(vg) \mathcal{M})=\frac{1}{\varepsilon}\{\mathcal{L}g-[v\partial_{x}(\rho\mathcal{ M}) \tag{2.22}\] \[-\partial_{x}(v(\rho\mathcal{M})).\mathcal{M}]-\rho(\partial_{t} \mathcal{M})\},\\ \partial_{t}\rho+\partial_{x}(v(\rho\mathcal{M}))+\varepsilon\partial_{x}(vg) =0,\\ -\partial_{x}^{2}\phi=\rho-h,\\ \langle g\rangle=0.\end{cases} \tag{2.22}\]

The equivalence between this model and the original VPFP system (2.1a)-(2.1b) has been demonstrated in [9]. Note that the conservation condition \(\langle g\rangle=0\) is usually satisfied by classical mesh-based AP schemes, but it fails to hold by DNN-predicated methods, unless incorporated into the loss function, as done in [17].

Now, the formal derivation of the high-field limit model can also be done easily from the micro-macro model (2.22). Taking the limit as \(\varepsilon\to 0\) in the second equation of the micro-macro model (2.22) yields:

\[\partial_{t}\rho+\partial_{x}\langle v(\rho\mathcal{M})\rangle=0. \tag{2.23}\]

Substituting the explicit expression for \(\mathcal{M}\) in the equation (2.23) reformulates it as:

\[\partial_{t}\rho-\partial_{x}(\rho\partial_{x}\phi)=0. \tag{2.24}\]

This results in the limit model:

\[\begin{cases}\partial_{t}\rho-\partial_{x}(\rho\partial_{x}\phi)=0,\\ -\partial_{x}^{2}\phi=\rho-h.\end{cases} \tag{2.25}\]

It is exactly the high-field limit model (2.10) in one-dimension, which justifies the AP property of the micro-macro model.

Finally, as shown in [9, 25], the Initial and periodic boundary conditions imposed on \(f\) can be directly translated into boundary conditions for \(\rho\) and \(g\).
Two formulations of asymptotic-preserving neural networks (APNNs)

We seek to approximate the solution to the VPFP system by neural networks. Under the framework of PINNs and other neural network-based methods, loss functions encoding the PDE constraints are minimized to converge on solutions. PINNs is our starting point but vanilla DNN approximations face limitations in tackling multiscale problems where both microscopic and macroscopic dynamics are involved. To deal with this difficulty, asymptotic-preserving neural networks (APNNs) have been proposed [17]. The key idea is to design a loss with AP property. Specifically, the loss function automatically transits from kinetic system to it asymptotic limit system when scale parameter \(\varepsilon\to 0\). In practice, the original multiscale system is first reformulated into an equivalent model with the AP property. This modified form is then embedded into the vanilla PINN loss function, with initial and boundary conditions treated as regularization terms.

To improve multiscale VPFP approximations of vanilla PINNs, we propose two types of loss functions with AP property. The constructions of the novel AP loss functions are respectively inspired by the derivations of the high-field limit models in Sections 1.2 and 1.3.

Consider the VPFP system with initial and boundary conditions over a bounded domain as an approximation of the original system (2.1a)-(2.1b), in which \((t,x,v)\in\mathcal{T}\times\mathcal{D}\times\Omega\) and \(\Omega\) is symmetric in \(v\):

\[\begin{cases}\partial_{t}f+v\cdot\nabla_{x}f=\frac{1}{\varepsilon}\mathcal{L }f,&(t,x,v)\in\mathcal{T}\times\mathcal{D}\times\Omega,\\ -\Delta_{x}\phi(t,x)=\rho(t,x)-h(x),&(t,x)\in\mathcal{T}\times\mathcal{D},\\ \mathcal{B}f=F_{\rm B},&(t,x,v)\in\mathcal{T}\times\partial\mathcal{D}\times \Omega,\\ \mathcal{B}\phi=\Phi_{\rm B},&(t,x)\in\mathcal{T}\times\partial\mathcal{D}\\ \mathcal{I}f=f_{0},&(t,x,v)\in\{t=0\}\times\mathcal{D}\times\Omega,\\ \mathcal{I}\phi=\phi_{0},&(t,x)\in\{t=0\}\times\mathcal{D}\end{cases} \tag{3.1}\]

where \(F_{\rm B},f_{0},\Phi_{\rm B},\phi_{0}\) are given functions; \(\partial\mathcal{D}\) is the boundary of \(\mathcal{D}\); and \(\mathcal{B},\mathcal{I}\) are boundary and initial operators, respectively. In the following sections, the limitations of vanilla PINNs for the multiscale VPFP system will be explained. We will also provide the detailed constructions of the two proposed AP loss functions and formally prove their AP properties.

### The failure of PINNs to resolve small scales

In this section, the deficiency of vanilla PINN are demonstrated in the case with small scale parameter \(\varepsilon\). First we use DNNs to parameterize two functions \(f\), \(\phi\) and denote the set of neural network parameters by \(\eta\) for brevity. Then two neural networks are implemented:

\[f_{\gamma}^{\rm NN}(t,x):=\log\left(1+\exp\left(\vec{f}_{\gamma}^{\rm NN}(t,x )\right)\right)\approx f(t,x),\ \ \text{and}\ \ \phi_{\gamma}^{\rm NN}(t,x)\approx\phi(t,x). \tag{3.2}\]

\(\vec{f}_{\gamma}^{\rm NN}\) and \(\phi_{\gamma}^{\rm NN}\) are both fully-connected neural networks. Note that the setting for \(f_{\gamma}^{\rm NN}\) guarantees the non-negativity property of distribution function \(f\).

The vanilla PINN loss is constructed by the least square of the residuals of the original VPFP system,
combined with penalty terms encoding the initial and boundary conditions:

\[\begin{split}\mathcal{R}^{\varepsilon}_{\text{PINN}}& =\frac{\mu_{1}^{\text{Vlasov}}}{\left|\mathcal{T}\times\mathcal{D}\times \Omega\right|}\int_{\mathcal{T}}\int_{\mathcal{D}}\int_{\Omega}\left| \varepsilon\partial_{t}f_{\gamma}^{\text{NN}}+\varepsilon v\partial_{x}f_{ \gamma}^{\text{NN}}-\mathcal{L}f_{\gamma}^{\text{NN}}\right|^{2}\ \mathrm{d}v\ \mathrm{d}x\mathrm{d}t\\ &+\frac{\mu_{1}^{\text{Poisson}}}{\left|\mathcal{T}\times\mathcal{D }\right|}\int_{\mathcal{T}}\int_{\mathcal{D}}\left|-\partial_{x}^{2}\phi_{ \gamma}^{\text{NN}}-\left(\rho_{\gamma}^{\text{NN}}-h\right)\right|^{2}\ \mathrm{d}x\mathrm{d}t\\ &+\frac{\mu_{2}^{f}}{\left|\mathcal{T}\times\mathcal{D}\times \Omega\right|}\int_{\mathcal{T}}\int_{\partial\mathcal{D}}\int_{\Omega}\left| \mathcal{B}f_{\gamma}^{\text{NN}}-F_{\text{B}}\right|^{2}\ \mathrm{d}v\mathrm{d}x\mathrm{d}t\\ &+\frac{\mu_{2}^{\phi}}{\left|\mathcal{T}\times\mathcal{D}\Omega \right|}\int_{\mathcal{T}}\int_{\partial\mathcal{D}}\left|\mathcal{B}\phi_{ \gamma}^{\text{NN}}-\Phi_{\text{B}}\right|^{2}\mathrm{d}x\mathrm{d}t\\ &+\frac{\mu_{3}^{f}}{\left|\mathcal{D}\times\Omega\right|}\int_ {\mathcal{D}}\int_{\Omega}\left|\mathcal{I}f_{\gamma}^{\text{NN}}-f_{0}\right| ^{2}\ \mathrm{d}v\mathrm{d}x+\frac{\mu_{3}^{\phi}}{\left|\mathcal{D}\right|}\int_{ \mathcal{D}}\left|\mathcal{I}\phi_{\gamma}^{\text{NN}}-\phi_{0}\right|^{2} \mathrm{d}x,\end{split} \tag{3.3}\]

where \(\mu_{1}\), \(\mu_{2}\) and \(\mu_{3}\) are the penalty weights to be tuned. For brevity, a formal notation is introduced:

\[\mathcal{R}^{\varepsilon}_{\text{PINN}}:=\mu_{1}\mathcal{R}^{\varepsilon}_{ \text{residual}}+\mu_{2}\mathcal{R}^{\varepsilon}_{\text{bc}}+\mu_{3} \mathcal{R}^{\varepsilon}_{\text{ic}}. \tag{3.4}\]

We now examine whether this loss of PINN method possesses the AP property. Wee only need to focus on the first term of (3.3):

\[\begin{split}\mathcal{R}^{\varepsilon}_{\text{residual}}=& \frac{1}{\left|\mathcal{T}\times\mathcal{D}\times\Omega\right|} \int_{\mathcal{T}}\int_{\mathcal{D}}\int_{\Omega}\left|\varepsilon\partial_{t }f_{\gamma}^{\text{NN}}+\varepsilon v\partial_{x}f_{\gamma}^{\text{NN}}- \mathcal{L}f_{\gamma}^{\text{NN}}\right|^{2}\ \mathrm{d}v\ \mathrm{d}x\mathrm{d}t\\ &+\frac{1}{\left|\mathcal{T}\times\mathcal{D}\right|}\int_{ \mathcal{T}}\int_{\mathcal{D}}\left|-\partial_{x}^{2}\phi_{\gamma}^{\text{NN}} -\left(\rho_{\gamma}^{\text{NN}}-h\right)\right|^{2}\ \mathrm{d}x\mathrm{d}t.\end{split} \tag{3.5}\]

Taking \(\varepsilon\to 0\), formally this leads to:

\[\begin{split}\mathcal{R}^{\varepsilon}_{\text{residual}}=& \frac{1}{\left|\mathcal{T}\times\mathcal{D}\times\Omega\right|} \int_{\mathcal{T}}\int_{\mathcal{D}}\int_{\Omega}\left|\mathcal{L}f_{\gamma}^{ \text{NN}}\right|^{2}\ \mathrm{d}v\ \mathrm{d}x\mathrm{d}t\\ &+\frac{1}{\left|\mathcal{T}\times\mathcal{D}\right|}\int_{ \mathcal{T}}\int_{\mathcal{D}}\left|-\partial_{x}^{2}\phi_{\gamma}^{\text{NN}} -\left(\rho_{\gamma}^{\text{NN}}-h\right)\right|^{2}\ \mathrm{d}x\mathrm{d}t.\end{split} \tag{3.6}\]

It is the least square loss of the following system:

\[\begin{cases}\left[\mathcal{L}f=0,\right.\\ -\partial_{x}^{2}\phi=\rho-h.\end{cases} \tag{3.7}\]

At the small scaling of \(\varepsilon\), our deviation shows that the PINN loss is dominated by its leading order terms. We are actually solving the equation \(\mathcal{L}f=0\), which gives \(f=\rho\mathcal{M}\). Clearly, this equation fails to produce the desired high-field limit equation (2.25). Therefore, training neural networks with vanilla PINN loss may lead to an inaccurate solution when \(\varepsilon\) is small. This will be furthermore, demonstrated by numerical experiments in Section 3.

### The micro-macro decomposition based APNN method

We now present an APNN method based on the micro-macro decomposition for the VPFP system. The main idea is to utilize PINN to solve the micro-macro system (2.22) rather than the original system (2.1a)-(2.1b). DNNs are employed to parameterize functions \(\rho\), \(g\), \(\phi\). Three independent neural networks are adopted for three functions. The set of parameters and the parameterized functions are separately denoted
by \(\theta\) and \(\rho^{\rm NN}_{\theta},g^{\rm NN}_{\theta}\), \(\phi^{\rm NN}_{\theta}\). In detail, \(\rho(t,x)\) is parameterized as:

\[\rho^{\rm NN}_{\theta}(t,x):=\log\Big{(}1+\exp\Big{(}\tilde{\rho}^{\rm NN}_{ \theta}(t,x)\Big{)}\Big{)}\approx\rho(t,x), \tag{3.8}\]

which conserves the non-negative property of density \(\rho\). And \(\phi(t,x)\) is parameterized as:

\[\phi^{\rm NN}_{\theta}(t,x)\approx\phi(t,x). \tag{3.9}\]

Moreover, \(g(t,x,v)\) is parameterized as:

\[g^{\rm NN}_{\theta}(t,x,v)=\tilde{g}^{\rm NN}_{\theta}(t,x,v)-\frac{1}{| \Omega|}\int_{\Omega}\tilde{g}^{\rm NN}_{\theta}(t,x,v)dv\approx g(t,x,v). \tag{3.10}\]

Here \(\tilde{\rho}^{\rm NN}_{\theta}\) and \(\tilde{g}^{\rm NN}_{\theta}\) and \(\phi^{\rm NN}_{\theta}\) are all fully-connected neural networks.

The design of \(g^{\rm NN}_{\theta}\) in (3.10) automatically guarantees the conservation condition in the micro-macro model (2.22):

\[\left\langle g^{\rm NN}_{\theta}\right\rangle=0,\quad\forall t,x. \tag{3.11}\]

As demonstrated in [17], this representation of \(g^{\rm NN}_{\theta}\) improves the accuracy of predicted results, compared to the method of incorporating this conservation mechanism as a soft constraint in loss. To verify the conservation property of \(g^{\rm NN}_{\theta}\), first, \(\left\langle\tilde{g}^{\rm NN}_{\theta}\right\rangle:=\int_{\mathbb{R}}\tilde {g}^{\rm NN}_{\theta}dv\) is approximated by \(\int_{\Omega}\tilde{g}^{\rm NN}_{\theta}dv\) in bounded domain. Then one has:

\[\begin{split} g^{\rm NN}_{\theta}(t,x,v)&=\tilde {g}^{\rm NN}_{\theta}(t,x,v)-\frac{1}{|\Omega|}\int_{\Omega}\tilde{g}^{\rm NN }_{\theta}(t,x,v)dv\\ &\approx\tilde{g}^{\rm NN}_{\theta}(t,x,v)-\frac{1}{|\Omega|} \langle\tilde{g}^{\rm NN}_{\theta}(t,x,v)\rangle.\end{split} \tag{3.12}\]

Finally, the following can be obtained:

\[\left\langle g^{\rm NN}_{\theta}\right\rangle\approx\int_{\Omega}g^{\rm NN}_ {\theta}dv=\int_{\Omega}\tilde{g}^{\rm NN}_{\theta}dv-\frac{1}{|\Omega|}(\int _{\Omega}\tilde{g}^{\rm NN}_{\theta}dv)(\int_{\Omega}dv)=0 \tag{3.13}\]

which confirms the conservation property.

To construct the loss function for the APNN method, we take the least squares of the residuals from the micro-macro system (2.22):

\[\begin{split}\mathcal{R}^{e}_{\rm APNN}=&\frac{ \lambda^{\rm macro}_{1}}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int _{\mathcal{D}}\left|\partial_{t}\rho^{\rm NN}_{\theta}+\partial_{x}\langle v( \rho^{\rm NN}_{\theta}\mathcal{M}^{\rm NN}_{\theta})\rangle+e\partial_{x} \langle vg^{\rm NN}_{\theta}\rangle\right|^{2}\;\mathrm{d}x\mathrm{d}t\\ &+\frac{\lambda^{\rm macro}_{1}}{|\mathcal{T}\times\mathcal{D} \times\mathcal{D}|}\int_{\mathcal{T}}\int_{\mathcal{D}}\left|\varepsilon \partial_{t}g^{\rm NN}_{\theta}+\varepsilon v(\partial_{x}g^{\rm NN}_{\theta} -\partial_{x}\langle vg^{\rm NN}_{\theta}\rangle\mathcal{M}^{\rm NN}_{ \theta})\right.\\ &-[\mathcal{L}_{g^{\rm NN}_{\theta}}-v\partial_{x}\langle\rho^{\rm NN }_{\theta}\mathcal{M}^{\rm NN}_{\theta}\rangle+\partial_{x}\langle v(\rho^{\rm NN }_{\theta}\mathcal{M}^{\rm NN}_{\theta})\rangle\mathcal{M}^{\rm NN}_{\theta}- \rho^{\rm NN}_{\theta}(\partial_{t}\mathcal{M}^{\rm NN}_{\theta})]\right|^{2} \;\mathrm{d}v\mathrm{d}x\mathrm{d}t\\ &+\frac{\lambda^{\rm Poisson}_{1}}{|\mathcal{T}\times\mathcal{D}|} \int_{\mathcal{T}}\int_{\mathcal{D}}\left|-\partial_{x}^{2}\phi^{\rm NN}_{ \theta}-(\rho^{\rm NN}_{\theta}-h)\right|^{2}\;\mathrm{d}x\mathrm{d}t\\ &+\frac{\lambda^{f}_{2}}{|\mathcal{T}\times\partial\mathcal{D} \times\Omega|}\int_{\mathcal{T}}\int_{\partial\mathcal{D}}\int_{\Omega}\left| \mathcal{B}(\rho^{\rm NN}_{\theta}\mathcal{M}^{\rm NN}_{\theta}+\varepsilon g ^{\rm NN}_{\theta})-F_{\rm B}\right|^{2}\;\mathrm{d}v\mathrm{d}x\mathrm{d}t \\ &+\frac{\lambda^{\rm d}_{2}}{|\mathcal{T}\times\partial\mathcal{D} |}\int_{\mathcal{T}}\int_{\partial\mathcal{D}}\left|\mathcal{B}\phi^{\rm NN}_{ \theta}-\Phi_{\rm B}\right|^{2}\;\mathrm{d}x\mathrm{d}t\\ &+\frac{\lambda^{f}_{3}}{|\mathcal{D}\times\mathcal{D}|}\int_{ \mathcal{D}}\int_{\Omega}\left|\mathcal{I}(\rho^{\rm NN}_{\theta}\mathcal{M}^{ \rm NN}_{\theta}+\varepsilon g^{\rm NN}_{\theta})-f_{0}\right|^{2}\;\mathrm{ d}v\mathrm{d}x\\ &+\frac{\lambda^{\rm d}_{3}}{|\mathcal{D}|}\int_{\mathcal{D}} \int_{\Omega}\left|\mathcal{I}\phi^{\rm NN}_{\theta}-\phi_{0}\right|^{2}\; \mathrm{d}x,\end{split} \tag{3.14}\]
where \(\mathcal{M}_{\theta}^{\text{NN}}\) can be calculated from \(\rho_{\theta}^{\text{NN}}\), \(g_{\theta}^{\text{NN}}\) and \(\phi_{\theta}^{\text{NN}}\). The loss is formally denoted as:

\[\mathcal{R}_{\text{APNN}}^{e}:=\lambda_{1}\mathcal{R}_{\text{residual}}^{e}+ \lambda_{2}\mathcal{R}_{\text{bc}}^{e}+\lambda_{3}\mathcal{R}_{\text{ic}}^{e}. \tag{3.15}\]

Here, we consider the periodic BCs. As precised at the end of Section 1.3, ICs and the periodic BCs of \(\rho\) and \(g\) can be directly deduced from the IC and periodic BC imposed on \(f\). Hence the initial and boundary conditions for \(f\) in the bounded domain VPFP system (3.1) are reformulated as:

\[\begin{cases}\mathcal{B}\rho=P_{\text{B}},&(t,x)\in\mathcal{T}\times \partial\mathcal{D}\\ \mathcal{B}g=G_{\text{B}},&(t,x,v)\in\mathcal{T}\times\partial\mathcal{D}\times \Omega\\ \mathcal{B}\phi=\Phi_{\text{B}},&(t,x)\in\mathcal{T}\times\partial\mathcal{D} \\ T\rho=\rho_{0},&(t,x)\in\{t=0\}\times\mathcal{D}\\ Tg=g_{0},&(t,x,v)\in\{t=0\}\times\mathcal{D}\times\Omega\\ T\phi=\phi_{0},&(t,x)\in\{t=0\}\times\mathcal{D}\end{cases} \tag{3.16}\]

where \(P_{\text{B}},G_{\text{B}},\Phi_{\text{B}},\rho_{0},\phi_{0}\) can be obtained from the corresponding settings of \(f\). In this case, the loss (3.14) is rewritten as:

\[\begin{split}\mathcal{R}_{\text{APNN}}^{e}&=\lambda_ {1}\mathcal{R}_{\text{residual}}^{e}+\frac{\lambda_{2}^{e}}{|\mathcal{T} \times\partial\mathcal{D}|}\int_{\mathcal{T}}\int_{\partial\mathcal{D}}\left| \mathcal{B}\mathcal{P}_{\theta}^{\text{NN}}-P_{\text{B}}\right|^{2}\,\text{d} x\text{d}t\\ &+\frac{\lambda_{2}^{g}}{|\mathcal{T}\times\partial\mathcal{D} \times\Omega\,|}\int_{\mathcal{T}}\int_{\partial\mathcal{D}}\left|\mathcal{B} \mathcal{B}_{\theta}^{\text{NN}}-G_{\text{B}}\right|^{2}\,\text{d}\text{vd}x \text{d}t+\frac{\lambda_{2}^{\theta}}{|\mathcal{T}\times\partial\mathcal{D}|} \int_{\mathcal{T}}\int_{\partial\mathcal{D}}\left|\mathcal{B}\mathcal{\phi}_{ \theta}^{\text{NN}}-\Phi_{\text{B}}\right|^{2}\text{d}x\text{d}t\\ &+\frac{\lambda_{3}^{e}}{|\mathcal{D}|}\int_{\mathcal{D}}\int_{ \Omega}\left|\mathcal{I}\rho_{\theta}^{\text{NN}}-\rho_{0}\right|^{2}\text{d} x+\frac{\lambda_{3}^{g}}{|\mathcal{D}\times\Omega|}\int_{\mathcal{D}}\int_{ \Omega}\left|\mathcal{I}g_{\theta}^{\text{NN}}-g_{0}\right|^{2}\,\text{d} \text{vd}x\\ &+\frac{\lambda_{3}^{\theta}}{|\mathcal{D}|}\int_{\mathcal{D}}\int_{ \Omega}\left|\mathcal{I}\phi_{\theta}^{\text{NN}}-\phi_{0}\right|^{2}\text{d} x.\end{split} \tag{3.17}\]

Wee now formally show the AP property of the loss function (3.14). Consider small \(\varepsilon\), we only need to focus on the residual term:

\[\begin{split}\mathcal{R}_{\text{residual}}^{e}&= \frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_{\mathcal{D}} \left|\partial_{r}\rho_{\theta}^{\text{NN}}+\partial_{x}\langle v(\rho_{ \theta}^{\text{NN}}\mathcal{M}_{\theta}^{\text{NN}})\rangle+\varepsilon \partial_{x}\langle vg_{\theta}^{\text{NN}}\rangle\right|^{2}\,\text{d}x \text{d}t\\ &+\frac{1}{|\mathcal{T}\times\mathcal{D}\times\Omega|}\int_{ \mathcal{T}}\int_{\mathcal{D}}\int_{\mathcal{\Omega}}\left|\,\varepsilon \partial_{t}g_{\theta}^{\text{NN}}+\varepsilon(\vartheta_{x}g_{\theta}^{\text{ NN}}-\partial_{x}\langle vg_{\theta}^{\text{NN}}\rangle\mathcal{M}_{ \theta}^{\text{NN}})\right.\\ &-\left[\mathcal{L}g_{\theta}^{\text{NN}}-v\partial_{x}\langle \rho_{\theta}^{\text{NN}}\mathcal{M}_{\theta}^{\text{NN}}\rangle+\partial_{x} \langle v(\rho_{\theta}^{\text{NN}}\mathcal{M}_{\theta}^{\text{NN}})\rangle \mathcal{M}_{\theta}^{\text{NN}}-\rho_{\theta}^{\text{NN}}(\partial_{r} \mathcal{M}_{\theta}^{\text{NN}})\right]|^{2}\,\text{d}v\text{d}x\text{d}t\\ &+\frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_{ \mathcal{D}}\left|-\partial_{x}^{2}\phi_{\theta}^{\text{NN}}-(\rho_{\theta}^{ \text{NN}}-h)\right|^{2}\,\text{d}x\text{d}t.\end{split} \tag{3.18}\]

Taking \(\varepsilon\to 0\), it formally becomes:

\[\begin{split}\mathcal{R}_{\text{residual}}^{e}&= \frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_{\mathcal{D}} \left|\partial_{\mu}\rho_{\theta}^{\text{NN}}+\partial_{x}\langle v(\rho_{ \theta}^{\text{NN}}\mathcal{M}_{\theta}^{\text{NN}})\rangle\right|^{2}\,\text{d} x\text{d}t\\ &+\frac{1}{|\mathcal{T}\times\mathcal{D}\times\Omega|}\int_{ \mathcal{T}}\int_{\mathcal{D}}\int_{\Omega}\left|-[\mathcal{L}g_{\theta}^{\text{ NN}}-v\partial_{x}\langle\rho_{\theta}^{\text{NN}}\mathcal{M}_{\theta}^{\text{NN}} \rangle+\partial_{x}\langle v(\rho_{\theta}^{\text{NN}}\mathcal{M}_{\theta}^{ \text{NN}})\rangle\mathcal{M}_{\theta}^{\text{NN}}-\rho_{\theta}^{\text{NN}}( \partial_{r}\mathcal{M}_{\theta}^{\text{NN}})\rangle\right|^{2}\,\text{d}v \text{d}x\text{d}t\\ &+\frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_{ \mathcal{D}}\left|-\partial_{x}^{2}\phi_{\theta}^{\text{NN}}-(\rho_{\theta}^{ \text{NN}}-h)\right|^{2}\,\text{d}x\text{d}t.\end{split} \tag{3.19}\]
This equation coincides with the least square loss of the system:

\[\begin{cases}\mathcal{L}g=[v\partial_{x}(\rho\mathcal{M})-\partial_{x}\langle v( \rho\mathcal{M})\rangle.\mathcal{M}]-\rho(\partial_{t}.\mathcal{M}),\\ \partial_{t}\rho+\partial_{x}\langle v(\rho\mathcal{M})\rangle=0,\\ -\partial_{x}^{2}\phi=\rho-h.\end{cases} \tag{3.20}\]

Similar to Section 1.3, the second equation in (3.20) produces \(\partial_{t}\rho-\partial_{x}\left(\rho\partial_{x}\phi\right)=0\). Combining it with \(-\partial_{x}^{2}\phi=\rho-h\), the high-field limit system (2.25) is derived. The AP property of this proposed DNN method based on the micro-macro decomposition is demonstrated. Therefore, this approach enables accurate predictions for both kinetic and high-field regimes.

### The mass conservation based APNN method

In this section, we propose another APNN method. This method is built upon an equivalent model of VPFP system, which incorporates the inherent mass conversion law to the origin VPFP system. Our new APNN method does not depend upon the explicit expression of the local Maxwellian, and nor it requires the initial data to be a local Maxwellian, thereby significantly broadening its potential applications.

#### 3.3.1 Enforcing the mass conservation law

A distinguished feature of this model is to enforce the local mass conversion law in the loss function. Let's starts from the origin VPFP equation (2.11):

\[\partial_{t}f+v\cdot\nabla_{x}f =\frac{1}{\varepsilon}\nabla_{v}\cdot[(v+\nabla_{x}\phi)f+\nabla _{v}f]=:\frac{1}{\varepsilon}\mathcal{L}f,\] \[-\Delta_{x}\phi =\rho-h(x).\]

Due to the property of mass conservation of Fokker-Planck-type operator \(\mathcal{L}\), one can integrate the first equation over \(v\) to obtain the equation of mass conservation law:

\[\partial_{t}\rho+\nabla_{x}\cdot\langle vf\rangle=0. \tag{3.21}\]

Our new model is the Vlasov equation (2.1a) together with an extra condition: the mass conservation law

\[\begin{cases}\partial_{t}f+v\cdot\nabla_{x}f=\frac{1}{\varepsilon}\mathcal{L }f,\\ \partial_{t}\rho+\nabla_{x}\cdot\langle vf\rangle=0,\\ -\Delta_{x}\phi=\rho-h,\\ \rho=\langle f\rangle.\end{cases} \tag{3.22}\]

This system is equivalent to the original system (2.1a)-(2.1b), since equation (3.21) is just the integration of the Vlasov equation (2.11) which is automatically satisfied in the continuous model but not necessarily so for DNN approximations. In the new model, the concerned macroscopic physical quantity \(\rho\) is naturally incorporated by the conservation law without any subsidiary conditions. And we will see the high-field limit equation naturally follows from the new model as \(\varepsilon\to 0\).

The formal derivation of the limit equation starts from the new model (3.22). Multiplying the first equation of the model (3.22) by \(\varepsilon\) and yielding \(\varepsilon\to 0\), this immediately gives \(\mathcal{L}f=0\). Since \(f\in\mathcal{N}(\mathcal{L})=\{f=\rho\mathcal{M}\), where \(\rho:=\langle f\rangle\}\), it follows that \(f=\rho\mathcal{M}\). Substituting this form of \(f\) into the second equation of
the model (3.22), and repeating the derivations in Section 1.3, limit system degenerates to:

\[\left\{\begin{array}{l}\partial_{t}\rho-\nabla_{x}\cdot(\rho\nabla_{x}\phi)=0, \\ -\triangle_{x}\phi=\rho-h(x).\end{array}\right. \tag{3.23}\]

Furthermore, one takes care of the corresponding initial-boundary value problem for the new model (3.22). The initial and boundary conditions of \(f\) are directly extended to \(\rho\), since \(\rho\) is direct integration of \(f\) over the velocity space.

#### 3.3.2 A new loss function with AP property

In this part, a new APNN method is proposed to solve the VPFP system. The key approach is incorporating the mass conservation based model (3.22), which possess the AP property, into the loss of the vanilla PINNs. The starting initial-boundary problem over a bounded domain is:

\[\begin{array}{ll}\left(\partial_{t}f+v\cdot\nabla_{x}f=\frac{1}{e}\mathcal{ L}f,&(t,x,v)\in\mathcal{T}\times\mathcal{D}\times\Omega\\ \partial_{t}\rho+\nabla_{x}\cdot(v\prime f)=0,&(t,x)\in\mathcal{T}\times \partial\mathcal{D}\\ -\triangle_{x}\phi(t,x)=\rho(t,x)-h(x),&(t,x)\in\mathcal{T}\times\mathcal{D},\\ \mathcal{B}f=F_{\text{B}},&(t,x,v)\in\mathcal{T}\times\partial\mathcal{D} \times\Omega\\ \mathcal{B}\rho=P_{\text{B}},&(t,x)\in\mathcal{T}\times\partial\mathcal{D} \\ \mathcal{B}\phi=\Phi_{\text{B}},&(t,x)\in\mathcal{T}\times\partial\mathcal{D} \\ Tf=f_{0},&(t,x,v)\in\{t=0\}\times\mathcal{D}\times\Omega\\ I\rho=\rho_{0},&(t,x)\in\{t=0\}\times\mathcal{D}\\ I\phi=\phi_{0},&(t,x)\in\{t=0\}\times\mathcal{D}\end{array} \tag{3.24}\]

where \(F_{\text{B}},f_{0},P_{\text{B}},\rho_{0},\Phi_{\text{B}},\phi_{0}\) are given functions.

First, we parametrize the three functions \(\rho\), \(f\) and \(\phi\) using DNNs, denoted by \(\eta\) for brevity. More precisely, we set:

\[\begin{array}{l}\rho_{\eta}^{\text{NN}}(t,x):=\log(1+\exp\left(\tilde{\rho} _{\eta}^{\text{NN}}(t,x)\right))\approx\rho(t,x),\\ f_{\eta}^{\text{NN}}(t,x):=\log(1+\exp\left(\tilde{f}_{\eta}^{\text{NN}}(t,x) \right))\approx f(t,x),\end{array} \tag{3.25}\]

which preserve the non-negative property of \(\rho\) and \(f\). The function \(\phi\) is parameterized as:

\[\phi_{\eta}^{\text{NN}}(t,x)\approx\phi(t,x). \tag{3.26}\]

Here \(\tilde{\rho}_{\theta}^{\text{NN}}\), \(\tilde{f}_{\theta}^{\text{NN}}\) and \(\phi_{\theta}^{\text{NN}}\) are all fully-connected neural networks. Next, we introduce the least square of the
residual of the mass conservation based model (3.22) as a new APNN loss:

\[\mathcal{R}^{e}_{\rm APNN}= \frac{\kappa_{1}^{\rm macro}}{|\mathcal{T}\times\mathcal{D}|}\int_ {\mathcal{T}}\int_{\mathcal{D}}\left|\partial_{\mu}\rho_{\eta}^{\rm NN}+ \partial_{x}\langle vf_{\eta}^{\rm NN}\rangle\right|^{2}\;\mathrm{d}x\mathrm{d}t \tag{3.27}\] \[+\frac{\kappa_{1}^{\rm kinetic}}{|\mathcal{T}\times\mathcal{D} \times\Omega|}\int_{\mathcal{T}}\int_{\mathcal{D}}\int_{\Omega}\left| \varepsilon\partial_{t}f_{\eta}^{\rm NN}+\varepsilon v\partial_{x}f_{\eta}^{ \rm NN}-\mathcal{L}f_{\eta}^{\rm NN}\right|^{2}\;\mathrm{d}v\;\mathrm{d}x \mathrm{d}t\] \[+\frac{\kappa_{1}^{\rm Poisson}}{|\mathcal{T}\times\mathcal{D}|} \int_{\mathcal{T}}\int_{\mathcal{D}}\left|-\partial_{x}^{2}\phi_{\eta}^{\rm NN }-\left(\rho_{\eta}^{\rm NN}-h\right)\right|^{2}\;\mathrm{d}x\mathrm{d}t\] \[+\frac{\kappa_{2}^{\rho}}{|\mathcal{T}\times\partial\mathcal{D}|} \int_{\mathcal{T}}\int_{\partial\mathcal{D}}\left|\mathcal{B}\rho_{\eta}^{\rm NN }-P_{\rm B}\right|^{2}\mathrm{d}x\mathrm{d}t\] \[+\frac{\kappa_{2}^{f}}{|\mathcal{T}\times\partial\mathcal{D}\times \Omega|}\int_{\mathcal{T}}\int_{\partial\mathcal{D}}\int_{\Omega}\left| \mathcal{B}f_{\eta}^{\rm NN}-F_{\rm B}\right|^{2}\;\mathrm{d}v\mathrm{d}x \mathrm{d}t\] \[+\frac{\kappa_{2}^{\rho}}{|\mathcal{T}\times\partial\mathcal{D}|} \int_{\mathcal{T}}\int_{\partial\mathcal{D}}\left|\mathcal{B}\phi_{\eta}^{\rm NN }-\Phi_{\rm B}\right|^{2}\mathrm{d}x\mathrm{d}t\] \[+\frac{\kappa_{3}^{\rho}}{|\mathcal{D}|}\int_{\mathcal{D}}\left| \mathcal{I}\rho_{\eta}^{\rm NN}-\rho_{0}\right|^{2}\mathrm{d}x+\frac{\kappa_{ 3}^{f}}{|\mathcal{D}\times\Omega|}\int_{\mathcal{D}}\int_{\Omega}\left| \mathcal{I}f_{\eta}^{\rm NN}-f_{0}\right|^{2}\;\mathrm{d}v\mathrm{d}x\] \[+\frac{\kappa_{4}}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{ T}}\int_{\mathcal{D}}\left|\left(f_{\eta}^{\rm NN}\right)-\rho_{\eta}^{\rm NN }\right|^{2}\;\mathrm{d}x\mathrm{d}t\]

For ease of description, loss \(\mathcal{R}^{e}_{\rm APNN}\) is abbreviated as:

\[\mathcal{R}^{e}_{\rm APNN}:=\kappa_{1}\mathcal{R}^{e}_{\rm residual}+\kappa_ {2}\mathcal{R}^{e}_{\rm bc}+\kappa_{3}\mathcal{R}^{e}_{\rm ic}. \tag{3.28}\]

At last, we show the AP property of the loss (3.27). We just focus on the residual term in loss:

\[\mathcal{R}^{e}_{\rm residual}= \frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_{ \mathcal{D}}\left|\partial_{\mu}\rho_{\eta}^{\rm NN}+\partial_{x}\langle vf_{ \eta}^{\rm NN}\rangle\right|^{2}\;\mathrm{d}x\mathrm{d}t \tag{3.29}\] \[+\frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_ {\mathcal{D}}\int_{\Omega}\left|\varepsilon\partial_{t}f_{\eta}^{\rm NN}+ \varepsilon v\partial_{x}f_{\eta}^{\rm NN}-\mathcal{L}f_{\eta}^{\rm NN} \right|^{2}\;\mathrm{d}v\;\mathrm{d}x\mathrm{d}t\] \[+\frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_ {\mathcal{D}}\left|-\partial_{x}^{2}\phi_{\eta}^{\rm NN}-\left(\rho_{\eta}^{ \rm NN}-h\right)\right|^{2}\;\mathrm{d}x\mathrm{d}t.\]

Taking \(\varepsilon\to 0\), formally it leads to:

\[\mathcal{R}^{e}_{\rm residual}= \frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_ {\mathcal{D}}\left|\partial_{\mu}\rho_{\eta}^{\rm NN}+\partial_{x}\langle vf_{ \eta}^{\rm NN}\rangle\right|^{2}\;\mathrm{d}x\mathrm{d}t \tag{3.30}\] \[+\frac{1}{|\mathcal{T}\times\mathcal{D}|}\int_{\mathcal{T}}\int_ {\mathcal{D}}\left|-\partial_{x}^{2}\phi_{\eta}^{\rm NN}-\left(\rho_{\eta}^{ \rm NN}-h\right)\right|^{2}\;\mathrm{d}x\mathrm{d}t.\]

This is the least square residual of the system:

\[\begin{cases}\mathcal{L}f=0,\\ \partial_{t}\rho+\partial_{x}\langle vf\rangle=0,\\ -\partial_{x}^{2}\phi=\rho-h.\end{cases} \tag{3.31}\]

The first equation \(\mathcal{L}f=0\) implies \(f\in\mathcal{N}(\mathcal{L})\), inducing \(f=\rho\mathcal{M}\). Plugging \(f=\rho\mathcal{M}\) into the second
equation and combining with the Poisson equation, one arrives at the limit system (2.25). Therefore, the AP property of the new APNN loss has been concluded.

Broadly, the mass conservation based APNN method holds several advantageous characteristics. First, the macroscopic physics quantity of interest is judiciously incorporated by the mass conservation law, free from restrictive assumptions. Secondly, the explicit Maxwellian form is rendered an unnecessary component. A wider applicability is promised, since not all kinetic-type equations possess explicit Maxwellians. Moreover, Our new method is applicable to both equilibrium and non-equilibrium initial data.

**Remark 3.1**.: _For the mass conservation based APNN method, the mass conservation condition \(\rho=\langle f\rangle\) is imposed to the loss as a soft constraint. This is different from the micro-macro decomposition based APNN method, in which the equivalent conservation condition \(\langle g\rangle=0\) is built as a hard constraint in DNNs. As reported in [20], the constraint-type for the conservation mechanism significantly effects the training accuracy of the DNNs. We will examine these two processing methods for the conservation conditions in the numerical experiment section._

### Empirical loss functions

Since the loss functions for APNNs and PINNs are defined by integrals, in practice, one needs to define the empirical loss as a Monte Carlo approximation to the integral loss. This involves the random selection to a small number of sub-domains (known as batches) to estimate the high-dimensional integrals. For the velocity-dependent operators \(\langle\cdot\rangle\) and \(\Pi(\cdot)\), the Gauss-Legendre quadrature provides additional accuracy. In detail, given \(\left\{w_{i},v_{i}^{\prime}\right\}_{i=1}^{n}\) are the nodes and weights, which \(\left\{v_{i}^{\prime}\right\}_{i=1}^{n}\) are the roots of the Legendre polynomials of degree \(n\) and \(\left\{w_{i}\right\}_{i=1}^{n}\) are the corresponding weights. Then an integral can be approximated, for example, \(\int_{-1}^{1}f(v)\mathrm{d}v\) by the summation of linear combination of \(f\left(v_{i}^{\prime}\right):\sum_{i-1}^{n}w_{i}f\left(v_{i}^{\prime}\right)\).

We start with the formulation of the empirical losses for vanilla PINN (3.3) defined as:

\[\begin{split}\mathcal{R}_{\mathrm{APNN}}^{e}=&\frac{ \mu_{1}^{\mathrm{Vlasov}}}{L_{1}}\sum_{l=1}^{L_{1}}\left|\delta\partial_{t}f_{ \gamma}^{\mathrm{NN}}(t_{l}^{r},x_{l}^{r},v_{l}^{r})+\varepsilon v_{l}^{r} \partial_{s}f_{\gamma}^{\mathrm{NN}}(t_{l}^{r},x_{l}^{r},v_{l}^{r})-\mathcal{ L}f_{\gamma}^{\mathrm{NN}}(t_{l}^{r},x_{l}^{r},v_{l}^{r})\right|^{2}\\ &+\frac{\mu_{1}^{\mathrm{Poisson}}}{L_{1}}\sum_{l=1}^{L_{1}}\left| -\partial_{x}^{2}\phi_{\gamma}^{\mathrm{NN}}(t_{l}^{r},x_{l}^{r})-(\rho_{ \gamma}^{\mathrm{NN}}(t_{l}^{r},x_{l}^{r})-h(x_{l}^{r}))\right|^{2}\\ &+\frac{\mu_{2}^{f}}{L_{2}}\sum_{l=1}^{L_{2}}\left|\mathcal{B}f_{ \gamma}^{\mathrm{NN}}(t_{l}^{b},x_{l}^{b},v_{l}^{b})-F_{\mathrm{B}}(t_{l}^{b},x_{l}^{b},v_{l}^{b})\right|^{2}+\frac{\mu_{2}^{\phi}}{L_{2}}\sum_{l=1}^{L_{2} }\left|\mathcal{B}\phi_{\gamma}^{\mathrm{NN}}(t_{l}^{b},x_{l}^{b})-\Phi_{ \mathrm{B}}(t_{l}^{b},x_{l}^{b})\right|^{2}\\ &+\frac{\mu_{3}^{f}}{L_{3}}\sum_{l=1}^{L_{3}}\left|\mathcal{I}f_{ \gamma}^{\mathrm{NN}}(t_{l}^{k},x_{l}^{k},v_{l}^{k})-f_{0}(t_{l}^{k},x_{l}^{k},v_{l}^{k})\right|^{2}+\frac{\mu_{3}^{\phi}}{L_{3}}\sum_{l=1}^{L_{1}}\left| \mathcal{I}\phi_{\gamma}^{\mathrm{NN}}(t_{l}^{k},x_{l}^{k})-\phi_{0}(t_{l}^{k},x_{l}^{k})\right|^{2}\end{split} \tag{3.32}\]

where \(\left\{(t_{l}^{r},x_{l}^{r},v_{l}^{r})\right\}_{l=1}^{L_{1}}\in\mathcal{T} \times\mathcal{D}\times\Omega\) is a batch of sample data for the residual, \(\left\{(t_{l}^{b},x_{l}^{b},v_{l}^{b})\right\}_{l=1}^{L_{2}}\in\mathcal{T} \times\partial\mathcal{D}\times\Omega\) is a batch of sample data, and \(\left\{(t_{l}^{k},x_{l}^{k},v_{l}^{k})\right\}_{l=1}^{L_{2}}\in\left\{0\right\} \times\mathcal{D}\times\Omega\) is a batch of sample data for initial conditions.

Similarly, the empirical loss function for the micro-macro decomposition based APNN method (3.14)
is defined as:

\[\begin{split}\mathcal{R}^{\text{ex}}_{\text{APNN}}=& \frac{\lambda_{1}^{\text{macro}}}{N_{1}}\sum_{i=1}^{N_{1}}\left| \partial_{t}\rho_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r})+\partial_{x}(v(\rho_ {\theta}^{\text{NN}}\mathcal{M}_{\theta}^{\text{NN}}))(t_{i}^{r},x_{i}^{r})+ \varepsilon\partial_{x}(vg_{\theta}^{\text{NN}})(t_{i}^{r},x_{i}^{r})\right|^ {2}\\ &+\frac{\lambda_{1}^{\text{micro}}}{N_{1}}\sum_{i=1}^{N_{1}}\left| \varepsilon\partial_{i}g_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r},v_{i}^{r})+ \varepsilon(v_{i}^{r}\partial_{x}g_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r},v _{i}^{r})-\partial_{x}(vg_{\theta}^{\text{NN}})(t_{i}^{r},x_{i}^{r})\mathcal{ M}_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r},v_{i}^{r})\right)\\ &-[\mathscr{L}_{\theta g}^{\text{NN}}(t_{i}^{r},x_{i}^{r},v_{i}^ {r})-v_{i}^{r}\partial_{x}(\rho_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r}) \mathcal{M}_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r},v_{i}^{r}))\\ &+\partial_{x}(v(\rho_{\theta}^{\text{NN}}\mathcal{M}_{\theta}^ {\text{NN}}))(t_{i}^{r},x_{i}^{r})\mathcal{M}_{\theta}^{\text{NN}}(t_{i}^{r},x _{i}^{r},v_{i}^{r})-\rho_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r})(\partial_{t }\mathcal{M}_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r},v_{i}^{r})))\right|^{2} \\ &\frac{\lambda_{1}^{\text{Poisson}}}{N_{1}}\sum_{i=1}^{N_{1}}\left|- \partial_{x}^{2}\phi_{\theta}^{\text{NN}}(t_{i}^{r},x_{i}^{r})-(\rho_{\theta} ^{\text{NN}}(t_{i}^{r},x_{i}^{r})-h(x_{i}^{r}))\right|^{2}\\ &+\frac{\lambda_{2}^{\prime}}{N_{2}}\sum_{i=1}^{N_{2}}\left| \mathcal{B}\left(\rho_{\theta}^{\text{NN}}(t_{i}^{b},x_{i}^{b})+\varepsilon g _{\theta}^{\text{NN}}(t_{i}^{b},x_{i}^{b},v_{i}^{b})\right)-F_{\text{B}}(t_{i} ^{b},x_{i}^{b},v_{i}^{b})\right|^{2}\\ &+\frac{\lambda_{2}^{\prime}}{N_{2}}\sum_{i=1}^{N_{2}}\left| \mathcal{B}\phi_{\theta}^{\text{NN}}(t_{i}^{b},x_{i}^{b})-\Phi_{\text{B}}(t_{i }^{b},x_{i}^{b})\right|^{2}\\ &+\frac{\lambda_{3}^{\prime}}{N_{3}}\sum_{i=1}^{N_{1}}\left| \mathcal{I}\left(\rho_{\theta}^{\text{NN}}\left(t_{i}^{b},x_{i}^{b}\right)+ \varepsilon g_{\theta}^{\text{NN}}\left(t_{i}^{b},x_{i}^{b},v_{i}^{b}\right) \right)-f_{0}\left(t_{i}^{b},x_{i}^{b},v_{i}^{b}\right)\right|^{2}\\ &+\frac{\lambda_{3}^{\prime}}{N_{3}}\sum_{i=1}^{N_{1}}\left| \mathcal{I}\phi_{\theta}^{\text{NN}}\left(t_{i}^{b},x_{i}^{b}\right)-\phi_{0} \left(t_{i}^{b},x_{i}^{b}\right)\right|^{2}\end{split} \tag{3.33}\]

where \(\{(t_{i}^{r},x_{i}^{r},v_{i}^{r})\}_{i=1}^{N_{1}}\in\mathcal{T}\times\mathcal{ D}\times\Omega\) is a batch of sample data for the residual, \(\{(t_{i}^{b},x_{i}^{b},v_{i}^{b})\}_{i=1}^{N_{2}}\in\mathcal{T}\times\partial \mathcal{D}\times\Omega\) is a batch of sample data for boundary conditions, and \(\{(t_{i}^{k},x_{i}^{k},v_{i}^{k})\}_{i=1}^{N_{1}}\in\{0\}\times\mathcal{D} \times\Omega\) is a batch of sample data for initial conditions.

The empirical loss function for the mass conservation based APNN method (3.27) is analogously de
fined as:

\[\begin{split}\mathcal{R}^{\text{e}}_{\text{APNN}}=& \frac{\kappa_{1}^{\text{macro}}}{M_{1}}\sum_{j=1}^{M_{1}}\left| \partial\rho_{\eta}^{\text{NN}}(t_{j}^{r},x_{j}^{r})+\partial_{x}(v_{\eta}^{ \text{NN}}(t_{j}^{r},x_{j}^{r}))\right|^{2}\\ &+\frac{\kappa_{1}^{\text{kinetic}}}{M_{1}}\sum_{j=1}^{M_{1}} \left|\varepsilon\partial_{t}f_{\eta}^{\text{NN}}(t_{j}^{r},x_{j}^{r},v_{j}^{ r})+\varepsilon v_{j}^{r}\partial_{x}f_{\eta}^{\text{NN}}(t_{j}^{r},x_{j}^{r},v_{j}^{ r})-\mathcal{L}f_{\eta}^{\text{NN}}(t_{j}^{r},x_{j}^{r},v_{j}^{r})\right|^{2}\\ &+\frac{\kappa_{1}^{\text{Poisson}}}{M_{1}}\sum_{j=1}^{M_{1}} \left|-\partial_{x}^{2}\phi_{\eta}^{\text{NN}}(t_{j}^{r},x_{j}^{r})-(\rho_{ \eta}^{\text{NN}}(t_{j}^{r},x_{j}^{r})-h(x_{j}^{r}))\right|^{2}\\ &+\frac{\kappa_{2}^{\rho}}{M_{2}}\sum_{j=1}^{M_{2}}\left| \mathcal{B}\rho_{\eta}^{\text{NN}}(t_{j}^{b},x_{j}^{b})-P_{\text{B}}(t_{j}^{b },x_{j}^{b})\right|^{2}+\frac{\kappa_{2}^{f}}{M_{2}}\sum_{j=1}^{M_{2}}\left| \mathcal{B}f_{\eta}^{\text{NN}}(t_{j}^{b},x_{j}^{b},v_{j}^{b})-F_{\text{B}}(t_ {j}^{b},x_{j}^{b},v_{j}^{b})\right|^{2}\\ &+\frac{\kappa_{2}^{\rho}}{M_{2}}\sum_{j=1}^{M_{2}}\left| \mathcal{B}\phi_{\eta}^{\text{NN}}(t_{j}^{b},x_{j}^{b})-\Phi_{\text{B}}(t_{j}^ {b},x_{j}^{b})\right|^{2}\\ &+\frac{\kappa_{3}^{\rho}}{M_{3}}\sum_{j=1}^{M_{3}}\left| \mathcal{I}\rho_{\eta}^{\text{NN}}(t_{j}^{k},x_{j}^{k})-\rho_{0}(t_{j}^{k},x_{ j}^{k})\right|^{2}+\frac{\kappa_{3}^{f}}{M_{3}}\sum_{j=1}^{M_{3}}\left| \mathcal{I}f_{\eta}^{\text{NN}}(t_{j}^{k},x_{j}^{k},v_{j}^{k})-f_{0}(t_{j}^{k},x_{j}^{k},v_{j}^{k})\right|^{2}\\ &+\frac{\kappa_{3}^{\phi}}{M_{3}}\sum_{j=1}^{M_{3}}\left| \mathcal{I}\phi_{\eta}^{\text{NN}}(t_{j}^{k},x_{j}^{k})-\phi_{0}(t_{j}^{k},x_{ j}^{k})\right|^{2}\\ &+\frac{\kappa_{4}}{M_{4}}\sum_{j=1}^{M_{4}}\left|\left\langle f_{ \eta}^{\text{NN}}\right\rangle(t_{j}^{r},x_{j}^{r})-\rho_{\eta}^{\text{NN}}(t_ {j}^{r},x_{j}^{r})\right|^{2}\end{split} \tag{3.34}\]

where \(\{(t_{j}^{r},x_{j}^{r},v_{j}^{r})\}_{j=1}^{M_{1}}\in\mathcal{T}\times\mathcal{ D}\times\Omega\) is a batch of sample data for the residual, \(\{(t_{j}^{b},x_{j}^{b},v_{j}^{b})\}_{j=1}^{M_{2}}\in\mathcal{T}\times\partial \mathcal{D}\times\Omega\) is a batch of sample data for boundary conditions, \(\{(t_{j}^{k},x_{j}^{k},v_{j}^{k})\}_{j=1}^{M_{1}}\in\{0\}\times\mathcal{D}\times\Omega\) is a batch of sample data for initial conditions, and \(\{(t_{j}^{r},x_{j}^{c})\}_{j=1}^{M_{4}}\in\mathcal{T}\times\mathcal{D}\) is a batch of sample data for the mass conservation.

With the empirical losses formed, one can generate the training set, and apply the optimization procedure for the empirical losses to approximate the solution of the VPFP system.

## 4 Numerical examples

In this section, extensive numerical results are presented for several problems chosen from kinetic regimes (\(\varepsilon\approx O(1)\)) to high-field regimes (\(\varepsilon\to 0\)), in order to verify two proposed APNN methods and compare their performances.

Some settings for the numerical experiments are stated as follows. The Adam version of the gradient descent algorithm is used for optimization of the loss functions. Neural network parameters \(\theta\), \(\eta\) and \(\gamma\) are initialized by the Xavier initialization. Fully-connected networks with 5 hidden layers and tanh activation function are used. The number of Gauss-quadrature points is set to be 32 for integration calculation. All hyper-parameters have been carefully tuned for optimal performance.

The reference solutions are obtained by standard finite difference methods. We will check the relative \(\ell^{2}\) errors of several macroscopic quantities, such as the density \(\rho(t,x)\) and the electric field \(E(t,x)=-\partial_{x}\phi(t,x)\), between the results of DNN methods and reference solutions. The root mean square errors (RMSE for brevity) are also provided when necessary. The specific definition of these two errors of \(\rho\) as an example
are:

\[\text{RMSE}(\rho):=\sqrt{\frac{1}{N}\sum_{j=1}^{N}\left|\rho_{\theta,j}^{\text{NN} }-\rho_{j}^{\text{ref}}\right|^{2}},\quad\ell^{2}(\rho):=\sqrt{\frac{\sum_{j} \left|\rho_{\theta,j}^{\text{NN}}-\rho_{j}^{\text{ref}}\right|^{2}}{\sum_{j} \left|\rho_{j}^{\text{ref}}\right|^{2}}}, \tag{4.1}\]

where \(N\) is the total number of mesh points.

### Problem I: Landau damping

In this part, the Landau damping case is exhibited to verify the AP property of the two proposed APNN loss functions from kinetic regimes to high-field regimes and compare their performances. The failure of vanilla PINN method in high-field regime is also demonstrated. The initial condition near the equilibrium takes the form:

\[f_{0}(x,v)=\frac{1}{\sqrt{2\pi}}\rho_{0}(x)\exp\left(-\frac{v^{2}}{2}\right), \quad(x,v)\in[0,2\pi/k]\times\mathbb{R}, \tag{4.2}\]

where the initial density \(\rho_{0}\) is:

\[\rho_{0}(x)=1+\alpha\cos(kx). \tag{4.3}\]

The wave number \(k\) is set to be \(0.5\) and the amplitude of the perturbation \(\alpha\) is set as \(0.05\). The computational domain is \([0,2\pi/k]\times[v_{\text{min}},v_{\text{max}}],-v_{\text{min}}=v_{\text{max }}=6\). The initial electric field \(\phi^{0}\) is the solution to the Poisson equation with \(h(x)=1\). It is solved analytically as:

\[\phi_{0}(x)=\frac{\alpha}{k^{2}}\cos(kx). \tag{4.4}\]

Periodic boundary conditions are applied in the \(x\) domain.

The detailed initial conditions (ICs) and boundary conditions (BCs) are specified for the APNN approaches. For the micro-macro decomposition based APNN method, the ICs for \(\rho_{0},g_{0},\phi_{0}\) are given by:

\[\rho_{0}(x)=1+\alpha\cos(kx),\quad\phi_{0}(x)=\frac{\alpha}{k^{2}}\cos(kx), \quad g_{0}(x,v)=\frac{1}{\varepsilon}\left\{f_{0}(x,v)-\rho_{0}(x)\mathcal{ M}_{0}(x,v)\right\}, \tag{4.5}\]

where \(\mathcal{M}_{0}\) is the Maxwellian as:

\[\mathcal{M}_{0}(x,v)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(v+\partial_{x} \phi_{0}(x))^{2}}{2}\right). \tag{4.6}\]

Periodic BCs on \(\rho,g\) are directly imposed based on the periodic BC for \(f\). For mass conservation based APNN method and vanilla PINN method, the IC and periodic BC on \(f\) are immediately translated into the related settings. To improve the numerical performance, exact periodic BCs are enforced. The ansatz is based on a Fourier basis. We also construct a transform \(\mathcal{T}:x\rightarrow\{\sin(kjx),\cos(kjx)\}_{j=1}^{n}\) before the first layer of neural networks and set \(n=1\).

Here, we present the performance of two APNN models in multiscale regimes, which are characterized by the parameter \(\varepsilon\). For each APNN model, we plot the density \(\rho\) and the electric field \(E\) as a function of \(x\) at different time, as well as the time evolution of the electric energy \(\|E(t)\|_{L^{2}}\). Figure 1 and Figure 2 illustrate that the two proposed APNN models achieve good agreement with the reference solutions across kinetic, intermediate, and high-field regimes. These numerical results confirm the AP property of two proposed APNN methods and supports our analysis. As a comparison, we plot the result of vanilla PINN method in Figure 3 and Figure 4. The vanilla PINN method exhibits slightly poorer performance in kinetic regimes (\(\varepsilon=1.0,0.5\)) and fails to capture the asymptotic limit regime. Overall, two proposed APNN models improve the vanilla PINN method and have the capability to reliably reproduce solutions for multiscale plasma physics systems.
Moreover, in the kinetic regime, we observe different performance of two proposed APNN methods at different time scales. The mass conservation based APNN method can capture the correct electric energy evolution over long duration up to \(t=5.0\) (see right column in Figure 2). However, the micro-macro decomposition based APNN method can only accurately characterize the dynamic behaviors in short duration up to \(t=1.0\) (see errors reported in Table 1, Table 2, Table 3). A possible reason is that the variable order of the non-equilibrium part \(g\) may influence the training of DNNs [32]. As time increases, the magnitude of \(g\) is significantly reduced. This may encumber the DNN training, thereby inducing the inaccuracy of the micro-macro decomposition based APNN method over long time duration. This weakness could be improved by employing self-adaptive PINNs [35, 42]. One promising strategy is to apply a time-dependent weight to non-equilibrium part \(g\)[35]. Another promising strategy is to allocate the time-dependent number of collocation points for non-equilibrium part \(g\)[42]. These self-adaptive strategies could assist in precisely characterizing \(g\) with variable order, thus improve the long time accuracy of the micro-macro decomposition based APNN method in VPFP system. In contrast, the mass conservation based APNN method gets rid of \(g\) and introduces the density \(\rho\). Density \(\rho\) and distribution \(f\) persist on equivalent first-orders, which broads its applicability over long time.

For the VPFP system, it is noteworthy that the hard or soft constraint-type of the mass conservation conditions do not significantly influence the accuracy of these two proposed APNN methods. In every numerical case tabled, the conservation condition (\(\langle g\rangle=0\)) is treated as a hard condition in the micro-macro decomposition based APNN method. We also examine the soft constraint (\(\langle f\rangle=\rho\)), constructed in the mass conservation based APNN method, which has a performance similar to its hard constraint counterpart (see Figure 1 and Figure 2). Thus, the APNN methods can achieve comparable performances with different constraint-type.
Figure 1: Multiscale linear Landau damping solved by the micro-macro decomposition based APNN method. Density \(\rho\) (left column) and electric field \(E\) (middle column) as functions of space \(x\) at \(t=0.0,0.5,1.0\); electric energy (right column) as a function of time \(t\) up to \(t=1.0\) with the kinetic and intermediate regimes, and \(t=5.0\) with the high-field regime. The electric energy is plotted in log scale. Neural networks are \([3,128,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,256,1]\) for \(g\). Batch size is 512 in domain and 256 for initial condition. Penalty \(\lambda_{1}=300\), \(\lambda_{3}=1\) for kinetic regime; penalty \(\lambda_{1}=0.5\), \(\lambda_{3}=1\) for intermediate regime; and penalty \(\lambda_{1}=0.5\), \(\lambda_{3}=1\) for high-field regime.
Figure 2: Multiscale linear Landau damping solved by the mass conservation based APNN method. Density \(\rho\) (left column) and electric field \(E\) (middle column) as functions of space \(x\) at \(t=0.0,0.5,1.0\); electric energy (right column) as a function of time \(t\) up to \(t=5.0\). Neural networks are \([3,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,256,1]\) for \(f\). Batch size is \(512\) in domain and \(256\) for initial condition and \(256\) for conservation condition \(\rho=\langle f\rangle\). Penalty \(\kappa_{1}=120\), \(\kappa_{3}=\kappa_{4}=1\) for kinetic regime; penalty \(\kappa_{1}=150\), \(\kappa_{3}=\kappa_{4}=1\) for intermediate regime; and penalty \(\kappa_{1}=500\), \(\kappa_{3}=\kappa_{4}=1\) for high-field regime.
\begin{table}

\end{table}
Table 2: Errors of linear Landau damping with intermediate regime (\(\varepsilon=0.5\)). Relative \(\ell^{2}\) error and RMSE of (a) density \(\rho\), (b) electric field \(E\) and (c) electric energy for the micro-macro decomposition based APNN method (MM for short), the mass conservation based APNN method (MC for short) and vanilla PINN method (PINN for short).

Figure 3: Multiscale linear Landau damping solved by vanilla PINN method with kinetic and intermediate regimes. Density \(\rho\) (column left) and electric field \(E\) (column middle) as functions of space \(x\) for \(t=0.0,0.5,1.0\); electric energy (right column) as a function of time \(t\) up to \(t=1.0\). Neural networks are \([3,128,128,128,128,128,128,1]\) for \(\phi\) and \([4,256,256,256,256,256,1]\) for \(f\). Batch size is 512 in domain and 256 for initial condition. Penalty \(\mu_{1}=30\), \(\mu_{3}=1\) for kinetic regime; and penalty \(\mu_{1}=50\), \(\mu_{3}=1\) for intermediate regime.
Figure 4: Linear Landau damping solved by vanilla PINN method with high-field regime (\(\varepsilon=0.01\)). Neural networks and batch size are the same as kinetic regime case. Penalty \(\mu_{1}=100\), \(\mu_{3}=1\).
### Problem II: The bump-on-tail case

In this part, we test whether proposed APNN methods are still effective when initial data is non-equilibrium. Consider the "double peak" function initial data:

\[f_{0}(x,v)=\frac{1}{\sqrt{2\pi}}\left(\frac{9}{10}\exp\left(-\frac{v^{2}}{2} \right)+\frac{2}{10}\exp\left(-4(v-4.5)^{2}\right)\right)(1+\alpha\cos(kx)),\]

where \(\alpha=0.05,k=0.5\) and \((x,v)\in[0,2\pi/k]\times[-v_{\max},v_{\max}]\) with \(v_{\max}=8\). The ICs for \(\rho,\phi\) are the same as the Landau damping case. The periodic BCs are applied.

Here, we demonstrate the performance of the two proposed models from kinetic regimes to high-field regimes. The micro-macro decomposition-based APNN method exhibits inaccuracy across all regimes, as depicted in Figure 5, Figure 6 and Figure 7. This inaccuracy may be attributable to the severe variations in the magnitude of \(g\), which occurs when the dynamics start from non-equilibrium initial data. As time progresses, these variation lead to a drastic change in order of \(g\), adversely impacting the performance of the neural network function \(g\). Consequently, this results in the inaccuracy of the learning outcomes [32]. Same as the micro-macro decomposition based APNN method for the Landau damping case, the inaccuracy of the learning outcomes attribute to the severe variations in the magnitude of \(g\). In contrast, the mass conservation based APNN method maintains the accurate outputs of neural networks for \(\rho\) and \(f\) at a consistent order of \(O(1)\), thereby ensuring good learning performance shown in Figure 8. As illustrated by this case, the APNN method exhibits adaptability to an non-equilibrium initial data, which underlines its flexibility and robustness.

\begin{table}

\end{table}
Table 3: Errors of linear Landau damping with high-field regime (\(\varepsilon=0.01\)). Relative \(\ell^{2}\) error and RMSE of (a) density \(\rho\), (b) electric field \(E\) and (c) electric energy for the micro-macro decomposition based APNN method (MM for short), the mass conservation based APNN method (MC for short) and vanilla PINN method (PINN for short).
Figure 5: The bump-on-tail case solved by the micro-macro decomposition based APNN method with kinetic regime (\(\varepsilon=1.0\)). Density \(\rho\) as a function of space \(x\) for \(t=0.0,0.5,1.0\) (top). Electric field \(E\) as a function of space \(x\) for \(t=0.0,0.5,1.0\) (bottom). Neural networks are \([3,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,256,1]\) for \(g\). Batch size is \(512\) in domain and \(256\) for initial condition. Penalty \(\lambda_{1}=50\) and \(\lambda_{3}=1\). Errors: Relative \(\ell^{2}\) error of \(\rho\) is \(1.93\times 10^{-3}\) at \(t=0.5\) and \(3.01\times 10^{-3}\) at \(t=1.0\). Relative \(\ell^{2}\) error of \(E\) is \(6.28\times 10^{-2}\) at \(t=0.5\) and \(1.58\times 10^{-1}\) at \(t=1.0\).

Figure 6: The bump-on-tail case solved by the micro-macro decomposition based APNN method with intermediate regime (\(\varepsilon=0.3\)). Density \(\rho\) as a function of space \(x\) for \(t=0.0,0.5,1.0\) (top). Electric field \(E\) as a function of space \(x\) for \(t=0.0,0.5,1.0\) (bottom). Neural networks and batch size are the same as kinetic regime case. Penalty \(\lambda_{1}=1\) and \(\lambda_{3}=1000\). Errors: Relative \(\ell^{2}\) error of \(\rho\) is \(4.76\times 10^{-3}\) at \(t=0.5\) and \(9.58\times 10^{-3}\) at \(t=1.0\). Relative \(\ell^{2}\) error of \(E\) is \(1.44\times 10^{-1}\) at \(t=0.5\) and \(7.86\times 10^{-1}\) at \(t=1.0\).
### Problem III: A Riemann problem

In this part, a one-dimensional Riemann problem is used to examine the performance of two proposed APNN methods in discontinuity cases. The piecewise constant initial condition is defined as:

\[\begin{cases}(\rho_{l},h_{l})=(1/8,1/2),&0\leq x<1/4,\\ \langle\rho_{m},h_{m}\rangle=(1/2,1/8),&1/4\leq x<3/4,\\ \langle\rho_{r},h_{r}\rangle=(1/8,1/2),&3/4\leq x\leq 1.\end{cases} \tag{4.7}\]

The exact initial condition of \(\phi\) is solved by the Poisson equation 2.1b as:

\[\begin{cases}\phi_{l}=3/16x^{2}-3/256,&0\leq x<1/4,\\ \phi_{m}=-3/16(x-1/2)^{2}+3/256,&1/4\leq x<3/4,\\ \phi_{r}=3/16(x-1)^{2}-3/256,&3/4\leq x\leq 1.\end{cases} \tag{4.8}\]

The corresponding first-order derivative \(\partial_{x}\phi\) is directly derived. Set \(f\) initially as:

\[f_{0}(x,v)=\frac{\rho_{0}}{\sqrt{2\pi}}\exp\left(-\frac{(v+\partial_{x}\phi_{0 }(x))^{2}}{2}\right). \tag{4.9}\]

Again periodic BC in \(x\) direction is applied.

We check the ability of two proposed APNN methods in tackling the discontinuities. Macroscopic variables \(\rho,\phi\) and flux \(j(t,x)=\int_{\mathcal{R}}vf(t,x,v)\mathrm{d}v\) are plotted in the high-field case (\(\varepsilon=0.001\)). Figure 9 and Figure 10 show that the proposed APNN methods can capture the sharp transition at singularities of the solution. High-frequency information such as sharp transitions are challenging to be learned by neural networks [46]. This difficulty can be covered by both APNN methods in this case.

We further note that the micro-macro decomposition based APNN method provides more accurate profiles in a shorter time frame, compared to the mass conservation based APNN method. This observation is supported by Table 4 and the bottom rows of Figure 9 and Figure 10. One possible explanation is due to the

Figure 7: The bump-on-tail case solved by the micro-macro decomposition based APNN method with high-field regime (\(\varepsilon=0.001\)). Density \(\rho\) as a function of space \(x\) for \(t=0.0,0.5,1.0\) (top). Electric field \(E\) as a function of space \(x\) for \(t=0.0,0.5,1.0\) (bottom). Neural networks, batch size and penalty are the same as intermediate regime case. Errors: Relative \(\ell^{2}\) error of \(\rho\) is \(2.94\times 10^{-3}\) at \(t=0.5\) and \(2.98\times 10^{-3}\) at \(t=1.0\). Relative \(\ell^{2}\) error of \(E\) is \(2.10\times 10^{-1}\) at \(t=0.5\) and \(3.63\times 10^{-1}\) at \(t=1.0\).
Figure 8: The bump-on-tail case solved by the mass conservation based APNN method. Density \(\rho\) (left column) and electric field \(E\) (right column) as functions of space \(x\) at \(t=0.0,0.5,1.0\). Neural networks are \([3,128,128,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,256,1]\) for \(f\). Batch size is \(512\) in domain, \(256\) for initial condition and \(256\) for conservation condition \(\rho=\langle f\rangle\). Penalty \(\kappa_{1}=50\), \(\kappa_{3}=\kappa_{4}=1\) for each regimes. Errors: \((a)\) Relative \(\ell^{2}\) error of \(\rho\) is \(1.51\times 10^{-3}\) at \(t=0.5\) and \(1.32\times 10^{-3}\) at \(t=1.0\). Relative \(\ell^{2}\) error of \(E\) is \(4.90\times 10^{-2}\) at \(t=0.5\) and \(6.53\times 10^{-2}\) at \(t=1.0\). \((b)\) Relative \(\ell^{2}\) error of \(\rho\) is \(1.02\times 10^{-3}\) at \(t=0.5\) and \(5.83\times 10^{-4}\) at \(t=1.0\). Relative \(\ell^{2}\) error of \(E\) is \(3.67\times 10^{-2}\) at \(t=0.5\) and \(3.22\times 10^{-2}\) at \(t=1.0\). \((c)\) Relative \(\ell^{2}\) error of \(\rho\) is \(1.93\times 10^{-3}\) at \(t=0.5\) and \(3.01\times 10^{-3}\) at \(t=1.0\). Relative \(\ell^{2}\) error of \(E\) is \(6.28\times 10^{-2}\) at \(t=0.5\) and \(1.58\times 10^{-1}\) at \(t=1.0\).
explicit representation of the non-equilibrium component \(g\) in the micro-macro decomposition based APNN method.This explicit representation enhances the DNNs' ability to accurately capture the non-equilibrium part of the distribution. The non-equilibrium part has a strong influence on the convergence of the governing VPFP system [32]. As a result, DNN could provide more detailed and insightful information. However, it's important to highlight that this improvement in accuracy is not universal. Specifically, it is only valid when the calculated time is relatively short and \(\rho\) and \(g\) are comparable in magnitude.

\begin{table}
\begin{tabular}{c c c c} \hline \hline \(\ell^{2}(\rho)\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(9.55\times 10^{-2}\) & \(7.60\times 10^{-2}\) \\ MC & \(7.36\times 10^{-2}\) & \(1.34\times 10^{-2}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline RMSE(\(\rho\)) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(3.46\times 10^{-2}\) & \(2.61\times 10^{-2}\) \\ MC & \(2.67\times 10^{-2}\) & \(4.59\times 10^{-2}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline \(\ell^{2}(E)\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.40\times 10^{-3}\) & \(1.58\times 10^{-3}\) \\ MC & \(8.32\times 10^{-4}\) & \(2.01\times 10^{-3}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.64\times 10^{-1}\) & \(1.46\times 10^{-1}\) \\ MC & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline RMSE(\(\text{flux}\)) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(3.23\times 10^{-3}\) & \(2.60\times 10^{-3}\) \\ MC & \(3.46\times 10^{-3}\) & \(4.17\times 10^{-3}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.64\times 10^{-1}\) & \(1.46\times 10^{-1}\) \\ MC & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.64\times 10^{-1}\) & \(1.46\times 10^{-1}\) \\ MC & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(3.23\times 10^{-3}\) & \(2.60\times 10^{-3}\) \\ MC & \(3.46\times 10^{-3}\) & \(4.17\times 10^{-3}\) \\ \hline \hline \end{tabular} \begin{tabular}{c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.0\) & \(t=0.2\) \\ \hline MM & \(1.76\times 10^{-1}\) & \(2.34\times 10^{-1}\) \\ \hline \hline \end{tabular} 
\begin{tabular}{c c} \hline \hline \(\ell^{2}(\text{flux})\) & \(t=0.
### Problem IV: Mixing regimes

The proposed APNN methods are now tested in the case with mixing regimes. In these regimes, the scale parameter \(\varepsilon\) fluctuates across several orders of magnitude in space. The parameter \(\varepsilon\) is constructed as:

\[\varepsilon(x)=\begin{cases}\varepsilon_{0}+\frac{1}{2}\left(\tanh(5-10x)+ \tanh(5+10x)\right),&-1\leq x\leq 0.3,\\ \varepsilon_{0},&0.3<x\leq 1,\end{cases} \tag{4.10}\]

where \(\varepsilon_{0}\) is set as \(0.001\). Hence both the kinetic and high-field regimes are under consideration. The initial conditions for \(\rho,f\) are given by:

\[\rho_{0}(x)=\frac{\sqrt{2\pi}}{6}(2+\sin(\pi x)),\quad f_{0}(x,v)=\frac{\rho_ {0}(x)}{\sqrt{2\pi}}\exp(-\frac{(v+\partial_{x}\phi_{0}(x))^{2}}{2}). \tag{4.11}\]

The initial electric field \(\phi_{0}(x)\) is analytically solved by the Poisson equation 2.1b with \(h(x)=\sqrt{2\pi}/3\). Periodic BC in \(x\) direction is imposed.

We report the capacity of the proposed APNN methods to resolve the case with mixing regimes. Figure 11 plots the scale parameter \(\varepsilon\) with a discontinuity at \(x=0.3\). The small \(\varepsilon\) induces transition of the solution is hard to be learned by DNNs-based methods as explained in the Riemann problem. Figure 12 and Table 5 show that both proposed APNN methods capture quite well the profile of density \(\rho\), except at \(x=0.3\) where the solution has a weak discontinuity. This inferior performance is ascribed to the inherent difficulty in learning high-frequency information by neural networks.

Figure 10: The Riemann problem by the mass conservation based APNN method in the high-field regime (\(\varepsilon=0.001\)). Density \(\rho\) as a function of space \(x\) at \(t=0.0,0.2\) (left column). Electric field \(E\) as a function of space \(x\) at \(t=0.0,0.2\) (middle column). Flux as a function of space \(x\) at \(t=0.0,0.2\) (right column). Neural networks are \([3,128,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,256,1]\) for \(f\). Batch size is \(512\) in domain, \(512\) for initial condition and \(256\) for conservation condition \(\rho=\langle f\rangle\). Penalty \(\kappa_{3}^{\rho}=3\), \(\kappa_{3}^{f}=1000\) and the rest of penalties equal to \(1\).
Figure 11: The scale parameter \(\varepsilon\) as a function of space \(x\) with several orders of magnitude.

Figure 12: The mixing regimes problem solved by the micro-macro decomposition based and mass conservation based APNN methods. Density \(\rho\) as a function of space \(x\) at \(t=0.0,0.1,0.2\). (\(a\)) Neural networks are \([3,128,128,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,1]\) for \(g\). Batch size is \(512\) in domain and \(256\) for initial condition. Penalty \(\lambda_{1}=0.1\) and \(\lambda_{3}=1\). (\(b\)) Neural networks are \([3,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,1]\) for \(f\). Batch size is \(512\) in domain, \(256\) for initial condition and \(256\) for conservation condition \(\rho=\langle f\rangle\). Penalty \(\kappa_{1}=\kappa_{3}=\kappa_{4}=1\).
### Problem V: The gravitational case

This section illustrates the efficacy of the proposed APNN methods for the VPFP system with gravitational force. Compared to the previous cases under repulsive electrostatic interaction, the attractive gravitational force induces different behavior. The existence of a unique weak solution locally in time has been proved in [36] for the limit equation with gravitational force.

The initial and boundary conditions applied are the same as in the Landau damping case. However, \(\phi_{0}\) is now the analytical solution to the equation \(\partial_{x}\phi_{0}(x)=\rho_{0}(x)-h(x)\) due to the gravitational force. As shown in Figure 13, the solutions obtained from both APNN methods exhibit good agreement with the exact solution. We conclude that both APNN methods correctly capture the dynamics of the VPFP system with gravitational force.

Figure 13: The gravitational case solved by the micro-macro decomposition based and mass conservation based APNN methods with high-field regime (\(\varepsilon=0.001\)). Density \(\rho\) (left column) and electric field \(E\) (right column) as functions of space \(x\) at \(t=0.0,0.2,0.5\). (\(a\)) Neural networks are \([3,128,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,256,1]\) for \(g\). Batch size is \(512\) in domain and \(256\) for initial condition. Penalty \(\lambda_{1}=\lambda_{3}=1\). (\(b\)) Neural networks are \([3,128,128,128,128,1]\) for \(\rho,\phi\) and \([4,256,256,256,256,256,1]\) for \(f\). Batch size is \(512\) in domain, \(256\) for initial condition and \(256\) for conservation condition \(\rho=\langle f\rangle\). Penalty \(\kappa_{1}=50\); and \(\kappa_{3}=\kappa_{4}=1\). Errors: (\(a\)) Relative \(\ell^{2}\) error of \(\rho\) is \(2.72\times 10^{-4}\) at \(t=0.2\) and \(9.78\times 10^{-4}\) at \(t=0.5\). Relative \(\ell^{2}\) error of \(E\) is \(3.00\times 10^{-3}\) at \(t=0.2\) and \(1.77\times 10^{-2}\) at \(t=0.5\). (\(b\)) Relative \(\ell^{2}\) error of \(\rho\) is \(5.82\times 10^{-4}\) at \(t=0.2\) and \(8.75\times 10^{-4}\) at \(t=0.5\). Relative \(\ell^{2}\) error of \(E\) is \(1.52\times 10^{-2}\) at \(t=0.2\) and \(1.94\times 10^{-2}\) at \(t=0.5\).

\begin{table}
\begin{tabular}{c c c c} \hline \hline \(\ell^{2}(\rho)\) & \(t=0.0\) & \(t=0.1\) & \(t=0.2\) \\ \hline MM & \(9.23\times 10^{-4}\) & \(1.36\times 10^{-2}\) & \(3.67\times 10^{-2}\) \\ MC & \(1.85\times 10^{-3}\) & \(2.06\times 10^{-2}\) & \(4.40\times 10^{-2}\) \\ \hline \hline \end{tabular}
\end{table}
Table 5: Error of the mixing regimes problem. Relative \(\ell^{2}\) error of \(\rho\) at \(t=0.0,0.1,0.2\) for the micro-macro decomposition based method (MM for short) and the mass conservation based method (MC for short).
### Problem VI: Uncertainty quantification (UQ) problems

This section highlights the potential of the proposed APNN methods for high-dimensional problems. Such problems present major challenges for classical mesh-based numerical schemes due to substantial computational costs. Here, we consider the uncertainty quantification (UQ) problems with high-dimensional random inputs. A typical setup of UQ is represented by the high-dimensional random input vectors. The initial conditions with random inputs are set as:

\[\rho_{0}(x,\mathbf{z})=\frac{\sqrt{2\pi}}{2}(2+\cos(2\pi x))+0.1\prod_{i=1}^{10 }\sin\left(\pi z_{i}\right),\ \ z=(z_{1},z_{2},\cdots,z_{10})\sim\mathcal{U}\left([-1,1]^{10}\right), \tag{4.12}\]

where \(z_{1},z_{2},\cdots,z_{10}\) are independent random variables following a uniform distribution over \([-1,1]\). The background charge \(h\) is set as:

\[h(x,\mathbf{z})=\sqrt{2\pi}+0.1\prod_{i=1}^{10}\sin\left(\pi z_{i}\right),\ \ z=(z_{1},z_{2},\cdots,z_{10})\sim\mathcal{U}\left([-1,1]^{10}\right). \tag{4.13}\]

Then one analytically solves initial \(\phi_{0}\) by the Poisson equation and defines:

\[f_{0}(x,v,\mathbf{z})=\frac{\rho_{0}}{\sqrt{2\pi}}\exp(-\frac{\left(v+\partial _{x}\phi_{0}\right)^{2}}{2}). \tag{4.14}\]

The computational domain is \((x,v)\in[0,1]\times[-6,6]\). Periodic boundary condition in \(x\) direction is applied. The high-field limit characterized by \(\varepsilon=0.001\) is considered.

We present the performance of the two proposed APNN methods in the high-dimensional case. Density \(\rho\) and electric field \(E\) are evaluated at \(t=0.0,0.05,0.1\) by taking expectation on \(10^{4}\) times simulations for \((z_{1},\cdots,z_{10})\). Figure 14 and Figure 15 exhibit good approximations between the APNN approximations and the exact solutions. These agreements imply the potential for APNN methods to resolve high-dimensional and multiscale problems.

## 5 Conclusion

In this work, we developed two numerical methods based on PINN-with tailored loss function suitable for solving the multiscale uncertain Vlasov-Poisson-Fokker-Planck system with possible the high-field regimes. Our work incorporates the Asymptotic-Preserving (AP) mechanism into the loss function. Besides the micro-macro decomposition based Asymptotic-Preserving neural network (APNN) method, we formulate another APNN method that enforces the mass conservation law. This method exhibits high adaptability for long time duration and non-equilibrium initial data. Diverse numerical examples are conducted across varied regimes to validate and compare the performance of both APNN methods. The numerical examples demonstrate the effectiveness of two proposed APNN methods for approximating multiscale VPFP systems. While the micro-macro decomposition based APNN method shows its accuracy in discountinuity cases, the mass conservation based method shows its broader applicability for the cases with long-duration or non-equilibrium initial data.

## Acknowledgement

Shi Jin is partially supported by the Strategic Priority Research Program of Chinese Academy of Sciences XDA25010401, the NSFC grant No. 12031013, the Shanghai Municipal Science and Technology Major
Figure 14: UQ problem solved by the micro-macro decomposition based APNN method. Density \(\rho\) (left column) and electric field \(E\) (right column) as functions of space \(x\) at \(t=0.0,0.05,0.1\). Neural networks are \([13,128,128,128,128,128,1]\) for \(\rho,\phi\) and \([14,256,256,256,256,256,1]\) for \(g\). Batch size is \(512\) in domain and \(256\) for initial condition. Penalty \(\lambda_{1}=50\), \(\lambda_{3}=1\) for kinetic regime; and \(\lambda_{1}=1\), \(\lambda_{3}=1\) for high-field regime. Errors: (\(a\)) Relative \(\ell^{2}\) error of \(\rho\) is \(4.71\times 10^{-3}\) at \(t=0.05\) and \(5.00\times 10^{-3}\) at \(t=0.1\). Relative \(\ell^{2}\) error of \(E\) is \(1.30\times 10^{-2}\) at \(t=0.05\) and \(2.28\times 10^{-2}\) at \(t=0.1\). (\(b\)) Relative \(\ell^{2}\) error of \(\rho\) is \(4.67\times 10^{-3}\) at \(t=0.05\) and \(7.86\times 10^{-3}\) at \(t=0.1\). Relative \(\ell^{2}\) error of \(E\) is \(2.11\times 10^{-2}\) at \(t=0.05\) and \(1.77\times 10^{-2}\) at \(t=0.1\).
Figure 15: UQ problem solved by the mass conservation based APNN method. Density \(\rho\) (left column) and electric field \(E\) (right column) as functions of space \(x\) at \(t=0.0,0.05,0.1\). Neural networks are \([13,128,128,128,128,128,1]\) for \(\rho,\phi\) and \([14,256,256,256,256,1]\) for \(f\). Batch size is \(512\) in domain, \(256\) for initial condition and \(256\) for conservation condition \(\rho=\langle f\rangle\). Penalty \(\kappa_{1}=100\), \(\kappa_{3}=\kappa_{4}=1\) for kinetic regime and \(\kappa_{1}=100\), \(\kappa_{3}=\kappa_{4}=1\) for high-field regime. Errors: \((a)\) Relative \(\ell^{2}\) error of \(\rho\) is \(7.83\times 10^{-3}\) at \(t=0.05\) and \(1.01\times 10^{-2}\) at \(t=0.1\). Relative \(\ell^{2}\) error of \(E\) is \(2.58\times 10^{-2}\) at \(t=0.05\) and \(4.02\times 10^{-2}\) at \(t=0.1\). \((b)\) Relative \(\ell^{2}\) error of \(\rho\) is \(7.57\times 10^{-3}\) at \(t=0.05\) and \(1.11\times 10^{-2}\) at \(t=0.1\). Relative \(\ell^{2}\) error of \(E\) is \(2.57\times 10^{-2}\) at \(t=0.05\) and \(3.34\times 10^{-2}\) at \(t=0.1\).
Project (2021SHZDZX0102), and the Fundamental Research Funds for the Central Universities. Zheng Ma is supported by NSFC Grant No. 12031013, No. 92270120 and Foundation of LCP.
