Introduction
The collisionality and magnetization of a plasma can vary by orders of magnitude within a single domain. Where the Knudsen number is large, the species require a kinetic description. Where it is small, they relax to local equilibrium and are governed by multi-fluid equations, which reduce in turn to ideal magnetohydrodynamics (MHD) in the limit of small normalized Larmor radius and Debye length [1 ] . Discrete-velocity methods that split transport and collisions must resolve the mean free path and the collision time to keep their numerical dissipation below the physical one [2 , 3 ] . Explicit particle-in-cell (PIC) methods must resolve the Debye length, the plasma period and the cyclotron period, and an explicit field solver also bounds the time step by the time in which light crosses a cell [4 ] . In the collision-dominated and MHD regimes, these scales are far smaller than those of the flow.
For the collisional scales, asymptotic-preserving (AP) schemes capture the continuum limit without resolving the collision time [3 , 5 ] . The unified gas-kinetic scheme (UGKS) [2 ] and the discrete unified gas-kinetic scheme (DUGKS) [6 ] construct their interface fluxes from the integral solution of the kinetic equation, which couples free transport and collisions within a time step. They thus capture the flow physics from free transport to the Navier–Stokes solutions according to the local cell Knudsen number, without resolving the mean free path and the collision time. Both have been extended to multicomponent plasmas coupled to Maxwell’s equations [7 , 8 ] . The unified gas-kinetic wave–particle (UGKWP) method [9 , 10 ] retains the multiscale flux of UGKS without its discrete velocity grid. Stochastic particles carry the part of the distribution function that streams freely over a time step, and an analytic hydrodynamic wave carries the part that collides within it. The particle fraction therefore decreases as the collision frequency increases and vanishes in the continuum limit.
For the electromagnetic scales, implicit PIC methods remove the stability limits that the plasma period and the light-crossing time impose on the time step [11 , 12 ] . The energy-conserving semi-implicit method (ECsim) [13 ] conserves the total energy exactly with a single linear system for the fields per time step, by building the current response from the same velocity update that then pushes the particles. Exact conservation of both the energy and the charge has been achieved in implicit electrostatic [14 ] and semi-implicit electromagnetic finite-element [15 ] PIC methods. As the Debye length tends to zero, the plasma becomes quasi-neutral, and Maxwell’s equations no longer determine the electric field. The reformulated Maxwell equation (RME) remains well posed in this limit, and particle methods based on it capture the quasi-neutral limit without resolving the plasma period and the Debye length [16 ] . AP particle methods have since been developed that conserve the discrete energy and satisfy Gauss’s law, either by correcting the electric field [17 ] or by conserving the charge [18 ] . For grid-based kinetic schemes, the RME has been used in a conservative semi-Lagrangian method that is AP in the quasi-neutral limit [19 ] and in DUGKS [8 ] for the Vlasov–Maxwell system with Bhatnagar–Gross–Krook (BGK) collisions [20 ] . The latter captures both the continuum and the quasi-neutral limits. With a reformulated Poisson equation, UGKS and UGKWP capture the same limits for electrostatic plasmas [21 ] .
The UGKWP method of Liu and Xu [1 ] advances Maxwell’s equations by an unconditionally stable Crank–Nicolson scheme and couples the momentum equations implicitly to the electromagnetic field. In the highly magnetized regime, its cell size and time step are then not restricted by the Debye length and the cyclotron period. The UGKWP method of Pu and Xu [22 ] advances Maxwell’s equations explicitly in their perfectly hyperbolic form [23 ] , with a time step limited by the speed of light, and couples the macroscopic variables of the species to the field by a Crank–Nicolson step. In UGKWP, the particles carry about the fraction e − Δ t / τ e^{-\Delta t/\tau} e − Δ t / τ Δ t \Delta t Δ t τ \tau τ τ \tau τ [1 ] and [22 ] , however, take the current from the macroscopic momentum equations and accelerate the particles in a separate update. Since the analytic wave carries part of the current, the velocity update that balances the exchange has to act on the particles and the wave alike.
Besides the energy balance, the discrete fields should satisfy Gauss’s law. With this coupling, the current in the field equation results from the velocity update, whereas the charge moves only in the transport steps. This current does not satisfy the discrete continuity equation, so that Gauss’s law does not hold by construction, unlike in PIC methods with charge-conserving current deposition [24 ] and in compatible particle–field discretizations [25 ] . Other methods control the Gauss defect by correcting the electric field [8 , 23 , 26 ] or the particle positions [27 ] . In quasi-neutral regions, a correction of the electric field at fixed charge can turn a small charge error into a large change of the field. Since part of the charge belongs to the wave, which is converted to and from particles in every transport step, a correction of the charge in UGKWP has to act on the particles and the wave together.
This paper presents UGKWP-EC, a semi-implicit UGKWP method for the multispecies BGK–Vlasov–Maxwell system, with an electromagnetic coupling that conserves energy. Between two UGKWP transport steps, a field step determines the electric field from an RME, a single linear system with a unique solution for any time step. Its current response is built from the implicit-midpoint velocity map, which then pushes the particles and the wave alike. Together with the implicit midpoint rule for Maxwell’s equations, this common map makes the field step conserve the discrete total energy exactly. At the end of each time step, a Gauss’s law correction adjusts the electric field and the masses of the particles and the wave while conserving the mass, momentum and energy of every charged species.
The field step remains stable without resolving the cyclotron period, the plasma period or the light-crossing time. The tests take the time step from a CFL condition on the flow and sound speeds of the species and use PIC, two-fluid and MHD references. In one dimension, they cover plasma waves, one of them with a collision frequency that varies over four orders of magnitude in space, and a Brio–Wu problem from kinetic plasmas to the MHD limit. In two dimensions, they cover the Orszag–Tang vortex and collisionless magnetic reconnection. In both, UGKWP-EC, whose time step does not depend on the speed of light, needs less CPU time than the UGKWP method of [22 ] . In the Orszag–Tang vortex, it also needs less CPU time than the implicit PIC code iPIC3D [12 ] with the same collision operator, and the saving grows as the collisions strengthen. In collisionless magnetic reconnection, it agrees with DUGKS [8 ] at less than one percent of the CPU time.
The remainder of the paper is organized as follows. § 2 states the BGK–Vlasov–Maxwell system, and § 3 reviews the UGKWP method. § 4 presents UGKWP-EC and its discrete properties, § 5 reports the numerical tests and the computational cost, and § 6 concludes the paper.
Governing equations
Consider N s N_s N s m s m_s m s q s q_s q s κ s = q s / m s \kappa_s=q_s/m_s κ s = q s / m s f s ( x , u , t ) f_s(\boldsymbol{x},\boldsymbol{u},t) f s ( x , u , t ) s s s x ∈ Ω ⊂ R d \boldsymbol{x}\in\Omega\subset\mathbb{R}^d x ∈ Ω ⊂ R d u ∈ R 3 \boldsymbol{u}\in\mathbb{R}^3 u ∈ R 3
∂ t f s + u ⋅ ∇ x f s + κ s ( E + u × B ) ⋅ ∇ u f s = M [ W s ] − f s τ s , s = 1 , … , N s , (1) \partial_t f_s+\boldsymbol{u}\cdot\nabla_{\boldsymbol{x}}f_s+\kappa_s(\boldsymbol{E}+\boldsymbol{u}\times\boldsymbol{B})\cdot\nabla_{\boldsymbol{u}}f_s=\frac{\mathcal{M}[\boldsymbol{W}_s]-f_s}{\tau_s},\qquad s=1,\ldots,N_s, \tag*{(1)} ∂ t f s + u ⋅ ∇ x f s + κ s ( E + u × B ) ⋅ ∇ u f s = τ s M [ W s ] − f s , s = 1 , … , N s , (1) where E \boldsymbol{E} E B \boldsymbol{B} B [20 ] , which relaxes f s f_s f s τ s \tau_s τ s M [ W s ] \mathcal{M}[\boldsymbol{W}_s] M [ W s ] s s s
W s = ∫ R 3 ψ f s d u = ( ρ s , P s , E s ) T , ψ ( u ) = ( 1 , u , 1 2 ∣ u ∣ 2 ) T . \boldsymbol{W}_s=\int_{\mathbb{R}^3}\psi f_s\,\mathrm{d}\boldsymbol{u}=(\rho_s,\boldsymbol{P}_s,\mathcal{E}_s)^{\mathsf{T}},\qquad\psi(\boldsymbol{u})=\left(1,\boldsymbol{u},\frac{1}{2}|\boldsymbol{u}|^2\right)^{\mathsf{T}}. W s = ∫ R 3 ψ f s d u = ( ρ s , P s , E s ) T , ψ ( u ) = ( 1 , u , 2 1 ∣ u ∣ 2 ) T . Its number density, macroscopic velocity and temperature are
n s = ρ s m s , U s = P s ρ s , T s = 2 3 n s k B ( E s − 1 2 ρ s ∣ U s ∣ 2 ) . n_s=\frac{\rho_s}{m_s},\qquad\boldsymbol{U}_s=\frac{\boldsymbol{P}_s}{\rho_s},\qquad T_s=\frac{2}{3n_sk_B}\left(\mathcal{E}_s-\frac{1}{2}\rho_s|\boldsymbol{U}_s|^2\right). n s = m s ρ s , U s = ρ s P s , T s = 3 n s k B 2 ( E s − 2 1 ρ s ∣ U s ∣ 2 ) . where k B k_B k B
M [ W s ] ( u ) = ρ s ( m s 2 π k B T s ) 3 / 2 exp ( − m s 2 k B T s ∣ u − U s ∣ 2 ) . \mathcal{M}[\boldsymbol{W}_s](\boldsymbol{u})=\rho_s\left(\frac{m_s}{2\pi k_BT_s}\right)^{3/2}\exp\left(-\frac{m_s}{2k_BT_s}|\boldsymbol{u}-\boldsymbol{U}_s|^2\right). M [ W s ] ( u ) = ρ s ( 2 π k B T s m s ) 3/2 exp ( − 2 k B T s m s ∣ u − U s ∣ 2 ) . The relaxation time of species s s s τ s = 1 / ν s \tau_s=1/\nu_s τ s = 1/ ν s ν s = χ s s n s \nu_s=\chi_{ss}n_s ν s = χ ss n s χ s s ≥ 0 \chi_{ss}\geq0 χ ss ≥ 0 s s s x \boldsymbol{x} x W s \boldsymbol{W}_s W s
The electromagnetic fields evolve according to Maxwell’s equations
ϵ 0 ∂ t E = μ 0 − 1 ∇ × B − J , ∂ t B = − ∇ × E , J = ∑ s κ s P s , (2) \epsilon_0\partial_t\boldsymbol{E}=\mu_0^{-1}\nabla\times\boldsymbol{B}-\boldsymbol{J},\qquad\partial_t\boldsymbol{B}=-\nabla\times\boldsymbol{E},\qquad\boldsymbol{J}=\sum_s\kappa_s\boldsymbol{P}_s, \tag*{(2)} ϵ 0 ∂ t E = μ 0 − 1 ∇ × B − J , ∂ t B = − ∇ × E , J = s ∑ κ s P s , (2) with the permittivity ϵ 0 \epsilon_0 ϵ 0 μ 0 \mu_0 μ 0 c = ( ϵ 0 μ 0 ) − 1 / 2 c=(\epsilon_0\mu_0)^{-1/2} c = ( ϵ 0 μ 0 ) − 1/2 ϵ 0 ∇ ⋅ E = ρ c \epsilon_0\nabla\cdot\boldsymbol{E}=\rho_c ϵ 0 ∇ ⋅ E = ρ c ρ c = ∑ s κ s ρ s \rho_c=\sum_s\kappa_s\rho_s ρ c = ∑ s κ s ρ s ∇ ⋅ B = 0 \nabla\cdot\boldsymbol{B}=0 ∇ ⋅ B = 0
∂ t ( ∑ s E s + ϵ 0 2 ∣ E ∣ 2 + 1 2 μ 0 ∣ B ∣ 2 ) + ∇ ⋅ ( ∑ s ∫ R 3 1 2 ∣ u ∣ 2 u f s d u + 1 μ 0 E × B ) = 0. \partial_t\left(\sum_s\mathcal{E}_s+\frac{\epsilon_0}{2}|\boldsymbol{E}|^2+\frac{1}{2\mu_0}|\boldsymbol{B}|^2\right)+\nabla\cdot\left(\sum_s\int_{\mathbb{R}^3}\frac{1}{2}|\boldsymbol{u}|^2\boldsymbol{u}f_s\,\mathrm{d}\boldsymbol{u}+\frac{1}{\mu_0}\boldsymbol{E}\times\boldsymbol{B}\right)=0. ∂ t ( s ∑ E s + 2 ϵ 0 ∣ E ∣ 2 + 2 μ 0 1 ∣ B ∣ 2 ) + ∇ ⋅ ( s ∑ ∫ R 3 2 1 ∣ u ∣ 2 u f s d u + μ 0 1 E × B ) = 0. The equations are nondimensionalized by the reference length x 0 x_0 x 0 t 0 t_0 t 0 ρ 0 \rho_0 ρ 0 m 0 m_0 m 0 B 0 B_0 B 0 u 0 = x 0 / t 0 u_0=x_0/t_0 u 0 = x 0 / t 0 E 0 = u 0 B 0 E_0=u_0B_0 E 0 = u 0 B 0
x = x 0 x ′ , t = t 0 t ′ , u = u 0 u ′ , f s = ρ 0 u 0 3 f s ′ , ( ρ s , P s , E s ) = ( ρ 0 ρ s ′ , ρ 0 u 0 P s ′ , ρ 0 u 0 2 E s ′ ) , \boldsymbol{x}=x_0\boldsymbol{x}',\qquad t=t_0t',\qquad\boldsymbol{u}=u_0\boldsymbol{u}',\qquad f_s=\frac{\rho_0}{u_0^3}f_s',\qquad(\rho_s,\boldsymbol{P}_s,\mathcal{E}_s)=(\rho_0\rho_s',\rho_0u_0\boldsymbol{P}_s',\rho_0u_0^2\mathcal{E}_s'), x = x 0 x ′ , t = t 0 t ′ , u = u 0 u ′ , f s = u 0 3 ρ 0 f s ′ , ( ρ s , P s , E s ) = ( ρ 0 ρ s ′ , ρ 0 u 0 P s ′ , ρ 0 u 0 2 E s ′ ) , m s = m 0 m s ′ , n s = ρ 0 m 0 n s ′ , T s = m 0 u 0 2 k B T s ′ , E = E 0 E ′ , B = B 0 B ′ . m_s=m_0m_s',\qquad n_s=\frac{\rho_0}{m_0}n_s',\qquad T_s=\frac{m_0u_0^2}{k_B}T_s',\qquad\boldsymbol{E}=E_0\boldsymbol{E}',\qquad\boldsymbol{B}=B_0\boldsymbol{B}'. m s = m 0 m s ′ , n s = m 0 ρ 0 n s ′ , T s = k B m 0 u 0 2 T s ′ , E = E 0 E ′ , B = B 0 B ′ . In these variables, (1), (2), Gauss’s law and ∇ ⋅ B = 0 \nabla\cdot\boldsymbol{B}=0 ∇ ⋅ B = 0 k B = 1 k_B=1 k B = 1
κ s ′ = κ s B 0 t 0 , ϵ 0 ′ = ϵ 0 B 0 2 ρ 0 , 1 μ 0 ′ = B 0 2 μ 0 ρ 0 u 0 2 , χ s s ′ = χ s s ρ 0 t 0 m 0 . \kappa_s'=\kappa_sB_0t_0,\qquad\epsilon_0'=\frac{\epsilon_0B_0^2}{\rho_0},\qquad\frac{1}{\mu_0'}=\frac{B_0^2}{\mu_0\rho_0u_0^2},\qquad\chi_{ss}'=\frac{\chi_{ss}\rho_0t_0}{m_0}. κ s ′ = κ s B 0 t 0 , ϵ 0 ′ = ρ 0 ϵ 0 B 0 2 , μ 0 ′ 1 = μ 0 ρ 0 u 0 2 B 0 2 , χ ss ′ = m 0 χ ss ρ 0 t 0 . The unified gas-kinetic wave–particle method
This section briefly reviews the UGKWP method [9 , 10 ] for a single species without electromagnetic fields, governed by the BGK equation
∂ t f + u ⋅ ∇ x f = g − f τ , g = M [ W ] , W = ∫ R 3 ψ f d u , (3) \partial_t f+\mathbf{u}\cdot\nabla_x f=\frac{g-f}{\tau},\qquad g=\mathcal{M}[\mathbf{W}],\qquad\mathbf{W}=\int_{\mathbb{R}^3}\psi f\,\mathrm{d}\mathbf{u}, \tag*{(3)} ∂ t f + u ⋅ ∇ x f = τ g − f , g = M [ W ] , W = ∫ R 3 ψ f d u , (3) with the relaxation time τ \tau τ
Finite-volume framework and interface flux
The domain is divided into cells Ω i \Omega_i Ω i ∣ Ω i ∣ |\Omega_i| ∣ Ω i ∣ S i j S_{ij} S ij Ω i \Omega_i Ω i Ω j \Omega_j Ω j j ∈ N ( i ) j\in N(i) j ∈ N ( i ) ∣ S i j ∣ |S_{ij}| ∣ S ij ∣ x i j \mathbf{x}_{ij} x ij n i j \mathbf{n}_{ij} n ij Ω j \Omega_j Ω j Δ t = t n + 1 − t n \Delta t=t^{n+1}-t^n Δ t = t n + 1 − t n W i \mathbf{W}_i W i F i j = ∫ t n t n + 1 ∫ R 3 u ⋅ n i j ψ f ( x i j , u , t ) d u d t \mathbf{F}_{ij}=\int_{t^n}^{t^{n+1}}\int_{\mathbb{R}^3}\mathbf{u}\cdot\mathbf{n}_{ij}\,\psi f(\mathbf{x}_{ij},\mathbf{u},t)\,\mathrm{d}\mathbf{u}\,\mathrm{d}t F ij = ∫ t n t n + 1 ∫ R 3 u ⋅ n ij ψ f ( x ij , u , t ) d u d t [10 ] . In UGKS [2 ] , the distribution function at the interface is the integral solution of (3) along the characteristics,
f ( x i j , u , t ) = 1 τ ∫ t n t e − ( t − t ′ ) / τ g ( x i j − u ( t − t ′ ) , u , t ′ ) d t ′ + e − ( t − t n ) / τ f n ( x i j − u ( t − t n ) , u ) , (4) f(\mathbf{x}_{ij},\mathbf{u},t)=\frac{1}{\tau}\int_{t^n}^{t}e^{-(t-t')/\tau}g\left(\mathbf{x}_{ij}-\mathbf{u}(t-t'),\mathbf{u},t'\right)\,\mathrm{d}t'+e^{-(t-t^n)/\tau}f^n\left(\mathbf{x}_{ij}-\mathbf{u}(t-t^n),\mathbf{u}\right), \tag*{(4)} f ( x ij , u , t ) = τ 1 ∫ t n t e − ( t − t ′ ) / τ g ( x ij − u ( t − t ′ ) , u , t ′ ) d t ′ + e − ( t − t n ) / τ f n ( x ij − u ( t − t n ) , u ) , (4) where f n f^n f n t n t^n t n x i j \mathbf{x}_{ij} x ij t n t^n t n g = g 0 + ( x − x i j ) ⋅ g x + ( t − t n ) g t g=g_0+(\mathbf{x}-\mathbf{x}_{ij})\cdot\mathbf{g}_x+(t-t^n)g_t g = g 0 + ( x − x ij ) ⋅ g x + ( t − t n ) g t Ω i \Omega_i Ω i u ⋅ n i j > 0 \mathbf{u}\cdot\mathbf{n}_{ij}>0 u ⋅ n ij > 0 Ω j \Omega_j Ω j u ⋅ n i j < 0 \mathbf{u}\cdot\mathbf{n}_{ij}<0 u ⋅ n ij < 0 f n = f 0 + ( x − x i j ) ⋅ f x f^n=f_0+(\mathbf{x}-\mathbf{x}_{ij})\cdot\mathbf{f}_x f n = f 0 + ( x − x ij ) ⋅ f x F i j = F i j e q + F i j f r \mathbf{F}_{ij}=\mathbf{F}_{ij}^{\mathrm{eq}}+\mathbf{F}_{ij}^{\mathrm{fr}} F ij = F ij eq + F ij fr
F i j e q = ∫ R 3 u ⋅ n i j ψ ( q 1 g 0 + q 2 u ⋅ g x + q 3 g t ) d u , F i j f r = ∫ R 3 u ⋅ n i j ψ ( q 4 f 0 + q 5 u ⋅ f x ) d u , (5) \begin{aligned}
\mathbf{F}_{ij}^{\mathrm{eq}}&=\int_{\mathbb{R}^3}\mathbf{u}\cdot\mathbf{n}_{ij}\,\psi\left(q_1g_0+q_2\mathbf{u}\cdot\mathbf{g}_x+q_3g_t\right)\,\mathrm{d}\mathbf{u},\\
\mathbf{F}_{ij}^{\mathrm{fr}}&=\int_{\mathbb{R}^3}\mathbf{u}\cdot\mathbf{n}_{ij}\,\psi\left(q_4f_0+q_5\mathbf{u}\cdot\mathbf{f}_x\right)\,\mathrm{d}\mathbf{u}, \tag*{(5)}
\end{aligned} F ij eq F ij fr = ∫ R 3 u ⋅ n ij ψ ( q 1 g 0 + q 2 u ⋅ g x + q 3 g t ) d u , = ∫ R 3 u ⋅ n ij ψ ( q 4 f 0 + q 5 u ⋅ f x ) d u , (5) with the coefficients
q 1 = Δ t − τ ( 1 − e − Δ t / τ ) , q 4 = τ ( 1 − e − Δ t / τ ) , q_1=\Delta t-\tau\left(1-e^{-\Delta t/\tau}\right),\qquad q_4=\tau\left(1-e^{-\Delta t/\tau}\right), q 1 = Δ t − τ ( 1 − e − Δ t / τ ) , q 4 = τ ( 1 − e − Δ t / τ ) , q 2 = 2 τ 2 ( 1 − e − Δ t / τ ) − τ Δ t − τ Δ t e − Δ t / τ , q 5 = τ Δ t e − Δ t / τ − τ 2 ( 1 − e − Δ t / τ ) , q_2=2\tau^2\left(1-e^{-\Delta t/\tau}\right)-\tau\Delta t-\tau\Delta t e^{-\Delta t/\tau},\qquad q_5=\tau\Delta t e^{-\Delta t/\tau}-\tau^2\left(1-e^{-\Delta t/\tau}\right), q 2 = 2 τ 2 ( 1 − e − Δ t / τ ) − τ Δ t − τ Δ t e − Δ t / τ , q 5 = τ Δ t e − Δ t / τ − τ 2 ( 1 − e − Δ t / τ ) , q 3 = Δ t 2 2 − τ Δ t + τ 2 ( 1 − e − Δ t / τ ) . q_3=\frac{\Delta t^2}{2}-\tau\Delta t+\tau^2\left(1-e^{-\Delta t/\tau}\right). q 3 = 2 Δ t 2 − τ Δ t + τ 2 ( 1 − e − Δ t / τ ) . The equilibrium at the interface is formed from the reconstructed conservative variables,
g 0 = M [ W 0 ] , W 0 = ∫ u ⋅ n i j > 0 ψ g L d u + ∫ u ⋅ n i j < 0 ψ g R d u , g_0=\mathcal{M}[\mathbf{W}_0],\qquad\mathbf{W}_0=\int_{\mathbf{u}\cdot\mathbf{n}_{ij}>0}\psi g_L\,\mathrm{d}\mathbf{u}+\int_{\mathbf{u}\cdot\mathbf{n}_{ij}<0}\psi g_R\,\mathrm{d}\mathbf{u}, g 0 = M [ W 0 ] , W 0 = ∫ u ⋅ n ij > 0 ψ g L d u + ∫ u ⋅ n ij < 0 ψ g R d u , W x = ∫ u ⋅ n i j > 0 ψ ∇ g L d u + ∫ u ⋅ n i j < 0 ψ ∇ g R d u , \mathbf{W}_x=\int_{\mathbf{u}\cdot\mathbf{n}_{ij}>0}\psi\nabla g_L\,\mathrm{d}\mathbf{u}+\int_{\mathbf{u}\cdot\mathbf{n}_{ij}<0}\psi\nabla g_R\,\mathrm{d}\mathbf{u}, W x = ∫ u ⋅ n ij > 0 ψ ∇ g L d u + ∫ u ⋅ n ij < 0 ψ ∇ g R d u , where g L g_L g L g R g_R g R x i j \mathbf{x}_{ij} x ij Ω i \Omega_i Ω i Ω j \Omega_j Ω j W x \mathbf{W}_x W x W 0 \mathbf{W}_0 W 0 g x \mathbf{g}_x g x g t g_t g t W x \mathbf{W}_x W x
∫ R 3 ψ g x d u = W x , ∫ R 3 ψ g t d u = − ∫ R 3 ψ u ⋅ g x d u . \int_{\mathbb{R}^3}\psi\mathbf{g}_x\,\mathrm{d}\mathbf{u}=\mathbf{W}_x,\qquad\int_{\mathbb{R}^3}\psi g_t\,\mathrm{d}\mathbf{u}=-\int_{\mathbb{R}^3}\psi\mathbf{u}\cdot\mathbf{g}_x\,\mathrm{d}\mathbf{u}. ∫ R 3 ψ g x d u = W x , ∫ R 3 ψ g t d u = − ∫ R 3 ψ u ⋅ g x d u . Wave–particle decomposition
The distribution function is represented in part by particles P k ( m k , x k , u k ) P_k(m_k,\mathbf{x}_k,\mathbf{u}_k) P k ( m k , x k , u k ) m k m_k m k x k \mathbf{x}_k x k u k \mathbf{u}_k u k (4) , a particle streams freely for a time t t t e − t / τ e^{-t/\tau} e − t / τ
t f , k = min ( − τ ln η k , Δ t ) , (6) t_{f,k}=\min\left(-\tau\ln\eta_k,\Delta t\right), \tag*{(6)} t f , k = min ( − τ ln η k , Δ t ) , (6) with η k \eta_k η k ( 0 , 1 ) (0,1) ( 0 , 1 ) x k + u k t f , k \mathbf{x}_k+\mathbf{u}_kt_{f,k} x k + u k t f , k Ω i \Omega_i Ω i
W i f r , p = ∑ x k + u k t f , k ∈ Ω i m k ψ ( u k ) − ∑ x k ∈ Ω i m k ψ ( u k ) . (7) \mathbf{W}_i^{\mathrm{fr},p}=\sum_{\mathbf{x}_k+\mathbf{u}_kt_{f,k}\in\Omega_i}m_k\boldsymbol{\psi}(\mathbf{u}_k)-\sum_{\mathbf{x}_k\in\Omega_i}m_k\boldsymbol{\psi}(\mathbf{u}_k). \tag*{(7)} W i fr , p = x k + u k t f , k ∈ Ω i ∑ m k ψ ( u k ) − x k ∈ Ω i ∑ m k ψ ( u k ) . (7) Particles with t f , k = Δ t t_{f,k}=\Delta t t f , k = Δ t t f , k < Δ t t_{f,k}<\Delta t t f , k < Δ t
W i h = W i − W i p , W i p = 1 ∣ Ω i ∣ ∑ x k ∈ Ω i m k ψ ( u k ) . \mathbf{W}_i^{h}=\mathbf{W}_i-\mathbf{W}_i^{p},\qquad\mathbf{W}_i^{p}=\frac{1}{|\Omega_i|}\sum_{\mathbf{x}_k\in\Omega_i}m_k\boldsymbol{\psi}(\mathbf{u}_k). W i h = W i − W i p , W i p = ∣ Ω i ∣ 1 x k ∈ Ω i ∑ m k ψ ( u k ) . Of the hydrodynamic wave, only the part
W i h p = e − Δ t / τ W i h (8) \mathbf{W}_i^{hp}=e^{-\Delta t/\tau}\mathbf{W}_i^{h} \tag*{(8)} W i h p = e − Δ t / τ W i h (8) moves freely over the whole step. This part is sampled from the Maxwellian of the reconstructed hydrodynamic wave into particles whose free-transport time is Δ t \Delta t Δ t ( 1 − e − Δ t / τ ) W i h (1-e^{-\Delta t/\tau})\mathbf{W}_i^{h} ( 1 − e − Δ t / τ ) W i h (5) of the wave without the part that the sampled particles carry,
F i j f r , w a v e = ∫ R 3 u ⋅ n i j ψ [ ( q 4 − Δ t e − Δ t / τ ) g 0 h + ( q 5 + Δ t 2 2 e − Δ t / τ ) u ⋅ g x h ] d u , (9) \mathbf{F}_{ij}^{\mathrm{fr},\mathrm{wave}}=\int_{\mathbb{R}^{3}}\mathbf{u}\cdot\mathbf{n}_{ij}\boldsymbol{\psi}\left[\left(q_4-\Delta t e^{-\Delta t/\tau}\right)g_0^{h}+\left(q_5+\frac{\Delta t^2}{2}e^{-\Delta t/\tau}\right)\mathbf{u}\cdot\mathbf{g}_x^{h}\right]\mathrm{d}\mathbf{u}, \tag*{(9)} F ij fr , wave = ∫ R 3 u ⋅ n ij ψ [ ( q 4 − Δ t e − Δ t / τ ) g 0 h + ( q 5 + 2 Δ t 2 e − Δ t / τ ) u ⋅ g x h ] d u , (9) where g 0 h g_0^{h} g 0 h g x h \mathbf{g}_x^{h} g x h x i j \mathbf{x}_{ij} x ij
W i n + 1 = W i n − 1 ∣ Ω i ∣ ∑ j ∈ N ( i ) ( F i j e q + F i j f r , w a v e ) ∣ S i j ∣ + W i f r , p ∣ Ω i ∣ , (10) \mathbf{W}_i^{n+1}=\mathbf{W}_i^n-\frac{1}{|\Omega_i|}\sum_{j\in N(i)}\left(\mathbf{F}_{ij}^{\mathrm{eq}}+\mathbf{F}_{ij}^{\mathrm{fr},\mathrm{wave}}\right)|S_{ij}|+\frac{\mathbf{W}_i^{\mathrm{fr},p}}{|\Omega_i|}, \tag*{(10)} W i n + 1 = W i n − ∣ Ω i ∣ 1 j ∈ N ( i ) ∑ ( F ij eq + F ij fr , wave ) ∣ S ij ∣ + ∣ Ω i ∣ W i fr , p , (10) where W i f r , p \mathbf{W}_i^{\mathrm{fr},p} W i fr , p
One time step consists of the following steps (Figure 1 ).
Figure 1. Composition of the distribution function in one step of the UGKWP method. The hatched region is the hydrodynamic wave, hollow circles are collisionless particles and solid circles are collisional particles. (a) Beginning of the step. (b) Classification of the particles by their free-transport times. (c) Sampling of W h p W^{hp} W h p 1 − e − Δ t / τ 1-e^{-\Delta t/\tau} 1 − e − Δ t / τ e − Δ t / τ e^{-\Delta t/\tau} e − Δ t / τ
Step 1. Form the hydrodynamic wave W i h = W i n − W i p \mathbf{W}_i^{h}=\mathbf{W}_i^n-\mathbf{W}_i^p W i h = W i n − W i p Figure 1 (a)).
Step 2. Obtain the free-transport time (6) of every remaining particle (Figure 1 (b)). Sample the collisionless particles of each cell from its hydrodynamic wave, with the expected moments ∣ Ω i ∣ W i h p |\Omega_i|\mathbf{W}_i^{hp} ∣ Ω i ∣ W i h p (8) and the free-transport time Δ t \Delta t Δ t Figure 1 (c)).
Step 3. Stream all particles freely for their free-transport times, accumulate the net particle moments W i f r , p \mathbf{W}_i^{\mathrm{fr},p} W i fr , p (7) , and compute the fluxes F i j e q \mathbf{F}_{ij}^{\mathrm{eq}} F ij eq (5) and F i j f r , w a v e \mathbf{F}_{ij}^{\mathrm{fr},\mathrm{wave}} F ij fr , wave (9) .
Step 4. Update the cell averages W i n + 1 \mathbf{W}_i^{n+1} W i n + 1 (10) and remove the collisional particles (Figure 1 (d)).
Wave–particle method for the BGK–Vlasov–Maxwell system
Time splitting and discrete state
The system (1), (2) is split into a transport–relaxation operator and a field operator,
A : ∂ t f s + u ⋅ ∇ x f s = M [ W s ] − f s τ s , ∂ t E = ∂ t B = 0 , V : ∂ t f s + κ s ( E + u × B ) ⋅ ∇ u f s = 0 , coupled to (2) . \begin{aligned}
\mathcal{A}: \quad&\partial_t f_s+\boldsymbol{u}\cdot\nabla_x f_s=\frac{\mathcal{M}[\boldsymbol{W}_s]-f_s}{\tau_s}, \qquad\partial_t\boldsymbol{E}=\partial_t\boldsymbol{B}=0,\\
\mathcal{V}: \quad&\partial_t f_s+\kappa_s(\boldsymbol{E}+\boldsymbol{u}\times\boldsymbol{B})\cdot\nabla_u f_s=0,\qquad\text{coupled to (2)}.
\end{aligned} A : V : ∂ t f s + u ⋅ ∇ x f s = τ s M [ W s ] − f s , ∂ t E = ∂ t B = 0 , ∂ t f s + κ s ( E + u × B ) ⋅ ∇ u f s = 0 , coupled to (2) . With δ = Δ t / 2 \delta=\Delta t/2 δ = Δ t /2 t n t^n t n t n + 1 = t n + Δ t t^{n+1}=t^n+\Delta t t n + 1 = t n + Δ t Figure 2 ): a transport step A δ \mathcal{A}_\delta A δ δ \delta δ § 4.2 ), the field step V Δ t \mathcal{V}_{\Delta t} V Δ t Δ t \Delta t Δ t § 4.3 ), a second transport step A δ \mathcal{A}_\delta A δ G \mathcal{G} G § 4.4 ). The stages after the first three substeps are labeled a a a b b b c c c
Figure 2. Substeps of one time step of UGKWP-EC, with the wave of a charged species at each stage. The field step maps the reconstructed wave h s a = M [ W s h , a ] h_s^a=\mathcal{M}[\boldsymbol{W}^{h,a}_s] h s a = M [ W s h , a ] h s b h_s^b h s b h s c = M [ W s h , c ] h_s^c=\mathcal{M}[\boldsymbol{W}^{h,c}_s] h s c = M [ W s h , c ] ( 1 + ζ s ) h s c (1+\zeta_s)h_s^c ( 1 + ζ s ) h s c h s n h_s^n h s n
At every stage, the discrete state of species s s s P k ( m k , x k , u k ) P_k(m_k,\boldsymbol{x}_k,\boldsymbol{u}_k) P k ( m k , x k , u k ) W s , i \boldsymbol{W}_{s,i} W s , i h s ( x , u ) h_s(\boldsymbol{x},\boldsymbol{u}) h s ( x , u ) W s , i h \boldsymbol{W}^{h}_{s,i} W s , i h E h \boldsymbol{E}_h E h B h \boldsymbol{B}_h B h § 4.3 ). The particles of species s s s Ω i \Omega_i Ω i P s , i P_{s,i} P s , i P s P_s P s W s , i p = ∣ Ω i ∣ − 1 ∑ k ∈ P s , i m k ψ ( u k ) \boldsymbol{W}^{p}_{s,i}=|\Omega_i|^{-1}\sum_{k\in P_{s,i}}m_k\boldsymbol{\psi}(\boldsymbol{u}_k) W s , i p = ∣ Ω i ∣ − 1 ∑ k ∈ P s , i m k ψ ( u k )
In the UGKWP method of § 3 , the hydrodynamic wave h s ( x , ⋅ ) = M [ W s h ( x ) ] h_s(\boldsymbol{x},\cdot)=\mathcal{M}[\boldsymbol{W}^{h}_s(\boldsymbol{x})] h s ( x , ⋅ ) = M [ W s h ( x )] W s h ( x ) \boldsymbol{W}^{h}_s(\boldsymbol{x}) W s h ( x ) W s , i h \boldsymbol{W}^{h}_{s,i} W s , i h h s a h_s^a h s a h s b h_s^b h s b h s c h_s^c h s c h s n + 1 h_s^{n+1} h s n + 1 W s , i h \boldsymbol{W}^{h}_{s,i} W s , i h
Transport step
The transport step A δ A_\delta A δ δ \delta δ (3) , with the relaxation time τ s \tau_s τ s n n n n + 1 n+1 n + 1 τ s , i \tau_{s,i} τ s , i W s , i n \mathbf{W}_{s,i}^n W s , i n
The fluxes of species s s s S i j S_{ij} S ij q 1 q_1 q 1 q 5 q_5 q 5 [28 ] , with the relaxation time τ s , i \tau_{s,i} τ s , i Ω i \Omega_i Ω i u n = u ⋅ n i j > 0 u_n=\mathbf{u}\cdot\mathbf{n}_{ij}>0 u n = u ⋅ n ij > 0 τ s , j \tau_{s,j} τ s , j Ω j \Omega_j Ω j u n < 0 u_n<0 u n < 0
F s , i j e q = ∫ u n > 0 u n ψ [ q 1 ( τ s , i ) g s , 0 + q 2 ( τ s , i ) u ⋅ ∇ g s , 0 + q 3 ( τ s , i ) g s , t ] d u + ∫ u n < 0 u n ψ [ q 1 ( τ s , j ) g s , 0 + q 2 ( τ s , j ) u ⋅ ∇ g s , 0 + q 3 ( τ s , j ) g s , t ] d u . \mathbf{F}_{s,ij}^{\mathrm{eq}} =
\int_{u_n>0} u_n\boldsymbol{\psi}\left[q_1(\tau_{s,i})g_{s,0}+q_2(\tau_{s,i})\mathbf{u}\cdot\nabla g_{s,0}+q_3(\tau_{s,i})g_{s,t}\right]\mathrm{d}\mathbf{u}
+\int_{u_n<0} u_n\boldsymbol{\psi}\left[q_1(\tau_{s,j})g_{s,0}+q_2(\tau_{s,j})\mathbf{u}\cdot\nabla g_{s,0}+q_3(\tau_{s,j})g_{s,t}\right]\mathrm{d}\mathbf{u}. F s , ij eq = ∫ u n > 0 u n ψ [ q 1 ( τ s , i ) g s , 0 + q 2 ( τ s , i ) u ⋅ ∇ g s , 0 + q 3 ( τ s , i ) g s , t ] d u + ∫ u n < 0 u n ψ [ q 1 ( τ s , j ) g s , 0 + q 2 ( τ s , j ) u ⋅ ∇ g s , 0 + q 3 ( τ s , j ) g s , t ] d u . where the interface equilibrium g s , 0 g_{s,0} g s , 0 ∇ g s , 0 \nabla g_{s,0} ∇ g s , 0 g s , t g_{s,t} g s , t W s \mathbf{W}_s W s (9) , the free-transport flux of the unsampled part of the wave is
F s , i j f r , w a v e = ∫ u n > 0 u n ψ [ ( q 4 ( τ s , i ) − δ p s , i ) h s , L + ( q 5 ( τ s , i ) + δ 2 2 p s , i ) u ⋅ ∇ h s , L ] d u + ∫ u n < 0 u n ψ [ ( q 4 ( τ s , j ) − δ p s , j ) h s , R + ( q 5 ( τ s , j ) + δ 2 2 p s , j ) u ⋅ ∇ h s , R ] d u . \mathbf{F}_{s,ij}^{\mathrm{fr,wave}} =
\int_{u_n>0} u_n\boldsymbol{\psi}\left[(q_4(\tau_{s,i})-\delta p_{s,i})h_{s,L}+\left(q_5(\tau_{s,i})+\frac{\delta^2}{2}p_{s,i}\right)\mathbf{u}\cdot\nabla h_{s,L}\right]\mathrm{d}\mathbf{u}
+\int_{u_n<0} u_n\boldsymbol{\psi}\left[(q_4(\tau_{s,j})-\delta p_{s,j})h_{s,R}+\left(q_5(\tau_{s,j})+\frac{\delta^2}{2}p_{s,j}\right)\mathbf{u}\cdot\nabla h_{s,R}\right]\mathrm{d}\mathbf{u}. F s , ij fr , wave = ∫ u n > 0 u n ψ [ ( q 4 ( τ s , i ) − δ p s , i ) h s , L + ( q 5 ( τ s , i ) + 2 δ 2 p s , i ) u ⋅ ∇ h s , L ] d u + ∫ u n < 0 u n ψ [ ( q 4 ( τ s , j ) − δ p s , j ) h s , R + ( q 5 ( τ s , j ) + 2 δ 2 p s , j ) u ⋅ ∇ h s , R ] d u . where h s , L h_{s,L} h s , L ∇ h s , L \nabla h_{s,L} ∇ h s , L h s , R h_{s,R} h s , R ∇ h s , R \nabla h_{s,R} ∇ h s , R h s n h_s^n h s n x i j x_{ij} x ij Ω i \Omega_i Ω i Ω j \Omega_j Ω j p s , i = e − δ / τ s , i p_{s,i}=e^{-\delta/\tau_{s,i}} p s , i = e − δ / τ s , i Ω i \Omega_i Ω i (8) .
As in (10) , the averages are updated in conservation form,
W s , i n + 1 = W s , i n − 1 ∣ Ω i ∣ ∑ j ∈ N ( i ) ( F s , i j e q + F s , i j f r , w a v e + F s , i j f r , p ) ∣ S i j ∣ , (11) \mathbf{W}_{s,i}^{n+1}=\mathbf{W}_{s,i}^n-\frac{1}{|\Omega_i|}\sum_{j\in N(i)}\left(\mathbf{F}_{s,ij}^{\mathrm{eq}}+\mathbf{F}_{s,ij}^{\mathrm{fr,wave}}+\mathbf{F}_{s,ij}^{\mathrm{fr,p}}\right)|S_{ij}|, \tag*{(11)} W s , i n + 1 = W s , i n − ∣ Ω i ∣ 1 j ∈ N ( i ) ∑ ( F s , ij eq + F s , ij fr , wave + F s , ij fr , p ) ∣ S ij ∣ , (11) where the particles of species s s s F s , i j f r , p ∣ S i j ∣ \mathbf{F}_{s,ij}^{\mathrm{fr,p}}|S_{ij}| F s , ij fr , p ∣ S ij ∣ S i j S_{ij} S ij Ω i \Omega_i Ω i Ω j \Omega_j Ω j − ∑ j ∈ N ( i ) F s , i j f r , p ∣ S i j ∣ -\sum_{j\in N(i)}\mathbf{F}_{s,ij}^{\mathrm{fr,p}}|S_{ij}| − ∑ j ∈ N ( i ) F s , ij fr , p ∣ S ij ∣ (7) .
At the end of the step, the mass of the wave is ρ s , i h , n + 1 = ρ s , i n + 1 − ρ s , i p , n + 1 \rho_{s,i}^{h,n+1}=\rho_{s,i}^{n+1}-\rho_{s,i}^{p,n+1} ρ s , i h , n + 1 = ρ s , i n + 1 − ρ s , i p , n + 1 W s , i n + 1 − W s , i p , n + 1 \mathbf{W}_{s,i}^{n+1}-\mathbf{W}_{s,i}^{p,n+1} W s , i n + 1 − W s , i p , n + 1 W s , i h , n = W s , i n − W s , i p , n \mathbf{W}_{s,i}^{h,n}=\mathbf{W}_{s,i}^n-\mathbf{W}_{s,i}^{p,n} W s , i h , n = W s , i n − W s , i p , n e − δ / τ s , i W s , i n e^{-\delta/\tau_{s,i}}\mathbf{W}_{s,i}^n e − δ / τ s , i W s , i n W s , i n + 1 − e − δ / τ s , i W s , i n = ( 1 − e − δ / τ s , i ) W s , i n + 1 \mathbf{W}_{s,i}^{n+1}-e^{-\delta/\tau_{s,i}}\mathbf{W}_{s,i}^n=(1-e^{-\delta/\tau_{s,i}})\mathbf{W}_{s,i}^{n+1} W s , i n + 1 − e − δ / τ s , i W s , i n = ( 1 − e − δ / τ s , i ) W s , i n + 1 W s , i n + 1 \mathbf{W}_{s,i}^{n+1} W s , i n + 1 W ^ s , i n + 1 = ( ρ s , i h , n + 1 / ρ s , i n + 1 ) W s , i n + 1 \widehat{\mathbf{W}}_{s,i}^{n+1}=(\rho_{s,i}^{h,n+1}/\rho_{s,i}^{n+1})\mathbf{W}_{s,i}^{n+1} W s , i n + 1 = ( ρ s , i h , n + 1 / ρ s , i n + 1 ) W s , i n + 1 P s , i n + 1 − P s , i p , n + 1 \mathbf{P}_{s,i}^{n+1}-\mathbf{P}_{s,i}^{p,n+1} P s , i n + 1 − P s , i p , n + 1 P ^ s , i n + 1 \widehat{\mathbf{P}}_{s,i}^{n+1} P s , i n + 1 W s , i \mathbf{W}_{s,i} W s , i P ^ s , i n + 1 \widehat{\mathbf{P}}_{s,i}^{n+1} P s , i n + 1 W ^ s , i n + 1 \widehat{\mathbf{W}}_{s,i}^{n+1} W s , i n + 1
P s , i h , n + 1 = ( 1 − ϖ s , i ) ( P s , i n + 1 − P s , i p , n + 1 ) + ϖ s , i P ^ s , i n + 1 , \mathbf{P}_{s,i}^{h,n+1}=(1-\varpi_{s,i})(\mathbf{P}_{s,i}^{n+1}-\mathbf{P}_{s,i}^{p,n+1})+\varpi_{s,i}\widehat{\mathbf{P}}_{s,i}^{n+1}, P s , i h , n + 1 = ( 1 − ϖ s , i ) ( P s , i n + 1 − P s , i p , n + 1 ) + ϖ s , i P s , i n + 1 , E s , i h , n + 1 = E ^ s , i n + 1 + 1 2 ρ s , i h , n + 1 ( ∣ P s , i h , n + 1 ∣ 2 − ∣ P ^ s , i n + 1 ∣ 2 ) . (12) \mathcal{E}_{s,i}^{h,n+1}=\widehat{\mathcal{E}}_{s,i}^{n+1}+\frac{1}{2\rho_{s,i}^{h,n+1}}\left(\lvert\mathbf{P}_{s,i}^{h,n+1}\rvert^2-\lvert\widehat{\mathbf{P}}_{s,i}^{n+1}\rvert^2\right). \tag*{(12)} E s , i h , n + 1 = E s , i n + 1 + 2 ρ s , i h , n + 1 1 ( ∣ P s , i h , n + 1 ∣ 2 − ∣ P s , i n + 1 ∣ 2 ) . (12) with ϖ s , i = N c / ( N c + N s , i h ) \varpi_{s,i}=N_c/(N_c+N^h_{s,i}) ϖ s , i = N c / ( N c + N s , i h ) N s , i h = ∣ Ω i ∣ ρ s , i h , n + 1 / m r e f , s N^h_{s,i}=|\Omega_i|\rho^{h,n+1}_{s,i}/m_{\mathrm{ref},s} N s , i h = ∣ Ω i ∣ ρ s , i h , n + 1 / m ref , s m r e f , s m_{\mathrm{ref},s} m ref , s N c > 0 N_c>0 N c > 0 N s , i h ≪ N c N^h_{s,i}\ll N_c N s , i h ≪ N c ϖ s , i \varpi_{s,i} ϖ s , i
Field step
The field step V Δ t \mathcal{V}_{\Delta t} V Δ t Δ t \Delta t Δ t x k a x_k^a x k a t n t^n t n h s a ( x , ⋅ ) = M [ W s h , a ( x ) ] h_s^a(x,\cdot)=\mathcal{M}[W_s^{h,a}(x)] h s a ( x , ⋅ ) = M [ W s h , a ( x )] W s , i h , a W_{s,i}^{h,a} W s , i h , a W s h , a = ( ρ s h , a , P s h , a , E s h , a ) T W_s^{h,a}=(\rho_s^{h,a},P_s^{h,a},\mathcal{E}_s^{h,a})^{\mathsf T} W s h , a = ( ρ s h , a , P s h , a , E s h , a ) T I i \mathcal{I}_i I i Ω i \Omega_i Ω i
The electric field is continuous and piecewise multilinear, determined by its values at the nodes of the mesh, and the magnetic field is piecewise constant, determined by its values in the cells. With the vectors e \boldsymbol{e} e b \boldsymbol{b} b E h ( x ) = N E ( x ) e E_h(x)=N_E(x)\boldsymbol{e} E h ( x ) = N E ( x ) e B h ( x ) = N B ( x ) b B_h(x)=N_B(x)\boldsymbol{b} B h ( x ) = N B ( x ) b N E N_E N E λ ℓ x ^ \lambda_\ell\hat{x} λ ℓ x ^ λ ℓ y ^ \lambda_\ell\hat{y} λ ℓ y ^ λ ℓ z ^ \lambda_\ell\hat{z} λ ℓ z ^ { λ ℓ } \{\lambda_\ell\} { λ ℓ } λ ℓ ≥ 0 \lambda_\ell\ge0 λ ℓ ≥ 0 ∑ ℓ λ ℓ = 1 \sum_\ell\lambda_\ell=1 ∑ ℓ λ ℓ = 1 Λ ( x ) \Lambda(x) Λ ( x ) λ ℓ ( x ) \lambda_\ell(x) λ ℓ ( x ) N B N_B N B x ^ \hat{x} x ^ y ^ \hat{y} y ^ z ^ \hat{z} z ^
M E = ∫ Ω N E T N E d x , M B = ∫ Ω N B T N B d x , C = ∫ Ω N B T ( ∇ × N E ) d x , M_E=\int_\Omega N_E^{\mathsf T}N_E\,\mathrm{d}x,\qquad M_B=\int_\Omega N_B^{\mathsf T}N_B\,\mathrm{d}x,\qquad C=\int_\Omega N_B^{\mathsf T}(\nabla\times N_E)\,\mathrm{d}x, M E = ∫ Ω N E T N E d x , M B = ∫ Ω N B T N B d x , C = ∫ Ω N B T ( ∇ × N E ) d x , the semi-discrete field equations with periodic or perfectly conducting boundaries are
ϵ 0 M E d e d t = μ 0 − 1 C T b − j , M B d b d t = − C e , j = ∫ Ω N E T J d x , \epsilon_0M_E\frac{\mathrm{d}\boldsymbol{e}}{\mathrm{d}t}=\mu_0^{-1}C^{\mathsf T}\boldsymbol{b}-\boldsymbol{j},\qquad M_B\frac{\mathrm{d}\boldsymbol{b}}{\mathrm{d}t}=-C\boldsymbol{e},\qquad\boldsymbol{j}=\int_\Omega N_E^{\mathsf T}J\,\mathrm{d}x, ϵ 0 M E d t d e = μ 0 − 1 C T b − j , M B d t d b = − C e , j = ∫ Ω N E T J d x , where j \boldsymbol{j} j
The electric and magnetic fields and the particle velocities are advanced by the implicit midpoint rule, with the rotation field B h r o t ( x ) = N B ( x ) b r o t B_h^{\mathrm{rot}}(x)=N_B(x)\boldsymbol{b}^{\mathrm{rot}} B h rot ( x ) = N B ( x ) b rot B B B B h r o t B_h^{\mathrm{rot}} B h rot [13 ] . For species s s s
L s ( x ) z = κ s z × B h r o t ( x ) , α s ( x ) = ( I − δ L s ( x ) ) − 1 , L_s(x)z=\kappa_s z\times B_h^{\mathrm{rot}}(x),\qquad\alpha_s(x)=\left(I-\delta L_s(x)\right)^{-1}, L s ( x ) z = κ s z × B h rot ( x ) , α s ( x ) = ( I − δ L s ( x ) ) − 1 , α s z = z + z × t s + t s ( t s ⋅ z ) 1 + ∣ t s ∣ 2 , t s = δ κ s B h r o t . \alpha_s z=\frac{z+z\times t_s+t_s(t_s\cdot z)}{1+|t_s|^2},\qquad t_s=\delta\kappa_s B_h^{\mathrm{rot}}. α s z = 1 + ∣ t s ∣ 2 z + z × t s + t s ( t s ⋅ z ) , t s = δ κ s B h rot . With e ˉ = 1 2 ( e n + e b ) \bar{\boldsymbol{e}}=\frac{1}{2}(\boldsymbol{e}^n+\boldsymbol{e}^b) e ˉ = 2 1 ( e n + e b ) k k k s s s
u k b − u k a = Δ t κ s ( N k e ˉ + u ˉ k × B k r o t ) , u ˉ k = 1 2 ( u k a + u k b ) , (13) u_k^b-u_k^a=\Delta t\ \kappa_s\left(N_k\bar{\boldsymbol{e}}+\bar{u}_k\times B_k^{\mathrm{rot}}\right),\qquad\bar{u}_k=\frac{1}{2}(u_k^a+u_k^b), \tag*{(13)} u k b − u k a = Δ t κ s ( N k e ˉ + u ˉ k × B k rot ) , u ˉ k = 2 1 ( u k a + u k b ) , (13) where N k = N E ( x k a ) N_k=N_E(x_k^a) N k = N E ( x k a ) B k r o t = B h r o t ( x k a ) B_k^{\mathrm{rot}}=B_h^{\mathrm{rot}}(x_k^a) B k rot = B h rot ( x k a )
u ˉ k = α k ( u k a + δ κ s N k e ˉ ) , u k b = 2 u ˉ k − u k a , (14) \bar{u}_k=\alpha_k\left(u_k^a+\delta\kappa_sN_k\bar{\boldsymbol{e}}\right),\qquad u_k^b=2\bar{u}_k-u_k^a, \tag*{(14)} u ˉ k = α k ( u k a + δ κ s N k e ˉ ) , u k b = 2 u ˉ k − u k a , (14) with α k = α s ( x k a ) \alpha_k = \alpha_s(x_k^a) α k = α s ( x k a ) x x x (14) assigns to a velocity u \mathbf{u} u h s a ( x , ⋅ ) h_s^a(x,\cdot) h s a ( x , ⋅ ) α s ( x ) ( u + δ κ s N E ( x ) e ˉ ) \alpha_s(x)(\mathbf{u}+\delta\kappa_sN_E(x)\bar{\mathbf{e}}) α s ( x ) ( u + δ κ s N E ( x ) e ˉ ) h s a ( x , ⋅ ) h_s^a(x,\cdot) h s a ( x , ⋅ )
P ‾ s h = α s ( P s h , a + δ κ s ρ s h , a N E e ‾ ) . \overline{\mathbf{P}}_s^h = \alpha_s\left(\mathbf{P}_s^{h,a}+\delta\kappa_s\rho_s^{h,a}N_E\overline{\mathbf{e}}\right). P s h = α s ( P s h , a + δ κ s ρ s h , a N E e ) . With the charge m k κ s m_k\kappa_s m k κ s k k k I i \mathcal{I}_i I i
j ‾ = ∑ s κ s ( ∑ k ∈ P s m k N k T u ‾ k + ∑ i I i [ N E T P ‾ s h ] ) . \overline{\mathbf{j}}=\sum_s\kappa_s\left(\sum_{k\in\mathcal{P}_s}m_kN_k^{\mathsf T}\overline{\mathbf{u}}_k+\sum_i\mathcal{I}_i\left[N_E^{\mathsf T}\overline{\mathbf{P}}_s^h\right]\right). j = s ∑ κ s ( k ∈ P s ∑ m k N k T u k + i ∑ I i [ N E T P s h ] ) . By (14) and the expression for P ‾ s h \overline{\mathbf{P}}_s^h P s h j ‾ \overline{\mathbf{j}} j e ‾ \overline{\mathbf{e}} e
j ^ p = ∑ s ∑ k ∈ P s m k κ s N k T α k u k a , X p = ∑ s ∑ k ∈ P s m k κ s 2 N k T α k N k , (15) \widehat{\mathbf{j}}_p=\sum_s\sum_{k\in\mathcal{P}_s}m_k\kappa_sN_k^{\mathsf T}\alpha_k\mathbf{u}_k^a,\qquad\mathbf{X}_p=\sum_s\sum_{k\in\mathcal{P}_s}m_k\kappa_s^2N_k^{\mathsf T}\alpha_kN_k, \tag*{(15)} j p = s ∑ k ∈ P s ∑ m k κ s N k T α k u k a , X p = s ∑ k ∈ P s ∑ m k κ s 2 N k T α k N k , (15) j ^ h = ∑ s ∑ i I i [ κ s N E T α s P s h , a ] , X h = ∑ s ∑ i I i [ κ s 2 ρ s h , a N E T α s N E ] , (16) \widehat{\mathbf{j}}_h=\sum_s\sum_i\mathcal{I}_i\left[\kappa_sN_E^{\mathsf T}\alpha_s\mathbf{P}_s^{h,a}\right],\qquad\mathbf{X}_h=\sum_s\sum_i\mathcal{I}_i\left[\kappa_s^2\rho_s^{h,a}N_E^{\mathsf T}\alpha_sN_E\right], \tag*{(16)} j h = s ∑ i ∑ I i [ κ s N E T α s P s h , a ] , X h = s ∑ i ∑ I i [ κ s 2 ρ s h , a N E T α s N E ] , (16) j ^ = j ^ p + j ^ h , X = X p + X h , j ‾ = j ^ + δ X e ‾ . \widehat{\mathbf{j}}=\widehat{\mathbf{j}}_p+\widehat{\mathbf{j}}_h,\qquad\mathbf{X}=\mathbf{X}_p+\mathbf{X}_h,\qquad\overline{\mathbf{j}}=\widehat{\mathbf{j}}+\delta\mathbf{X}\overline{\mathbf{e}}. j = j p + j h , X = X p + X h , j = j + δ X e . The field update is the implicit midpoint rule for the semi-discrete field equations, with the discrete current j ‾ \overline{\mathbf{j}} j
ϵ 0 M E e b − e n Δ t = μ 0 − 1 C T b ‾ − j ‾ , \epsilon_0M_E\frac{\mathbf{e}^b-\mathbf{e}^n}{\Delta t}=\mu_0^{-1}C^{\mathsf T}\overline{\mathbf{b}}-\overline{\mathbf{j}}, ϵ 0 M E Δ t e b − e n = μ 0 − 1 C T b − j , M B b b − b n Δ t = − C e ‾ , b ‾ = 1 2 ( b b + b n ) . (17) M_B\frac{\mathbf{b}^b-\mathbf{b}^n}{\Delta t}=-C\overline{\mathbf{e}},\qquad\overline{\mathbf{b}}=\frac{1}{2}\left(\mathbf{b}^b+\mathbf{b}^n\right). \tag*{(17)} M B Δ t b b − b n = − C e , b = 2 1 ( b b + b n ) . (17) The discrete Faraday law in (17) gives b ‾ = b n − δ M B − 1 C e ‾ \overline{\mathbf{b}}=\mathbf{b}^n-\delta M_B^{-1}C\overline{\mathbf{e}} b = b n − δ M B − 1 C e δ \delta δ e ‾ \overline{\mathbf{e}} e
( ϵ 0 M E + δ 2 X + δ 2 μ 0 C T M B − 1 C ) e ‾ = ϵ 0 M E e n − δ j ^ + δ μ 0 C T b n . (18) \left(\epsilon_0M_E+\delta^2\mathbf{X}+\frac{\delta^2}{\mu_0}C^{\mathsf T}M_B^{-1}C\right)\overline{\mathbf{e}}=\epsilon_0M_E\mathbf{e}^n-\delta\widehat{\mathbf{j}}+\frac{\delta}{\mu_0}C^{\mathsf T}\mathbf{b}^n. \tag*{(18)} ( ϵ 0 M E + δ 2 X + μ 0 δ 2 C T M B − 1 C ) e = ϵ 0 M E e n − δ j + μ 0 δ C T b n . (18) Equation (18) is an RME [16 , 19 ] , in which the current response (16) is the discrete form of the current equation ∂ t J = ∑ s κ s 2 ( ρ s E + P s × B ) \partial_t\mathbf{J}=\sum_s\kappa_s^2(\rho_s\mathbf{E}+\mathbf{P}_s\times\mathbf{B}) ∂ t J = ∑ s κ s 2 ( ρ s E + P s × B ) V \mathcal{V} V
The RME has exactly one solution for every Δ t > 0 \Delta t>0 Δ t > 0 Proposition 1 ), and the endpoint fields are
e b = 2 e ‾ − e n , b b = 2 b ‾ − b n = b n − Δ t M B − 1 C e ‾ . (19) \mathbf{e}^b=2\overline{\mathbf{e}}-\mathbf{e}^n,\qquad\mathbf{b}^b=2\overline{\mathbf{b}}-\mathbf{b}^n=\mathbf{b}^n-\Delta tM_B^{-1}C\overline{\mathbf{e}}. \tag*{(19)} e b = 2 e − e n , b b = 2 b − b n = b n − Δ t M B − 1 C e . (19) Since the field step changes b \mathbf{b} b − Δ t M B − 1 C e ‾ -\Delta tM_B^{-1}C\overline{\mathbf{e}} − Δ t M B − 1 C e (19) and the other substeps leave b \mathbf{b} b b n − b 0 \mathbf{b}^n-\mathbf{b}^0 b n − b 0 M B − 1 C M_B^{-1}C M B − 1 C ∫ Ω ∇ λ ℓ ⋅ B h d x \int_\Omega\nabla\lambda_\ell\cdot\mathbf{B}_h\,\mathrm{d}x ∫ Ω ∇ λ ℓ ⋅ B h d x ℓ \ell ℓ
With e ‾ \overline{\mathbf{e}} e (18) , every particle is pushed once by (14) . Written without u ‾ k \overline{\mathbf{u}}_k u k u k b = R s ( x k a ) u k a + d s ( x k a ) \mathbf{u}_k^b=R_s(x_k^a)\mathbf{u}_k^a+\mathbf{d}_s(x_k^a) u k b = R s ( x k a ) u k a + d s ( x k a )
R s ( x ) = 2 α s ( x ) − I , d s ( x ) = Δ t κ s α s ( x ) N E ( x ) e ‾ . (20) R_s(x)=2\alpha_s(x)-I,\qquad\mathbf{d}_s(x)=\Delta t\,\kappa_s\alpha_s(x)N_E(x)\overline{\mathbf{e}}. \tag*{(20)} R s ( x ) = 2 α s ( x ) − I , d s ( x ) = Δ t κ s α s ( x ) N E ( x ) e . (20) where R s = ( I + δ L s ) ( I − δ L s ) − 1 R_s=(I+\delta L_s)(I-\delta L_s)^{-1} R s = ( I + δ L s ) ( I − δ L s ) − 1 δ L s \delta L_s δ L s h s a ( x , ⋅ ) h_s^a(x,\cdot) h s a ( x , ⋅ ) u ↦ R s ( x ) u + d s ( x ) u\mapsto R_s(x)u+d_s(x) u ↦ R s ( x ) u + d s ( x )
h s b ( x , u ) = h s a ( x , R s ( x ) T [ u − d s ( x ) ] ) . (21) h_s^b(x,u)=h_s^a\left(x,R_s(x)^{\mathsf T}[u-d_s(x)]\right). \tag*{(21)} h s b ( x , u ) = h s a ( x , R s ( x ) T [ u − d s ( x )] ) . (21) The map moves the velocity U s h , a U_s^{h,a} U s h , a h s a ( x , ⋅ ) h_s^a(x,\cdot) h s a ( x , ⋅ ) R s U s h , a + d s R_sU_s^{h,a}+d_s R s U s h , a + d s R s R_s R s h s b ( x , ⋅ ) = M [ W s h , b ( x ) ] h_s^b(x,\cdot)=\mathcal{M}[W_s^{h,b}(x)] h s b ( x , ⋅ ) = M [ W s h , b ( x )]
ρ s h , b = ρ s h , a , U s h , b = R s U s h , a + d s , T s h , b = T s h , a , \rho_s^{h,b}=\rho_s^{h,a},\qquad U_s^{h,b}=R_sU_s^{h,a}+d_s,\qquad T_s^{h,b}=T_s^{h,a}, ρ s h , b = ρ s h , a , U s h , b = R s U s h , a + d s , T s h , b = T s h , a , P s h , b = R s P s h , a + ρ s h , a d s , E s h , b = E s h , a + d s ⋅ R s P s h , a + 1 2 ρ s h , a ∣ d s ∣ 2 . (22) P_s^{h,b}=R_sP_s^{h,a}+\rho_s^{h,a}d_s,\qquad\mathcal{E}_s^{h,b}=\mathcal{E}_s^{h,a}+d_s\cdot R_sP_s^{h,a}+\frac{1}{2}\rho_s^{h,a}|d_s|^2. \tag*{(22)} P s h , b = R s P s h , a + ρ s h , a d s , E s h , b = E s h , a + d s ⋅ R s P s h , a + 2 1 ρ s h , a ∣ d s ∣ 2 . (22) The cell averages are updated by the increments of the particles and of the wave,
Δ W s , i p = 1 ∣ Ω i ∣ ∑ k ∈ P s , i m k [ ψ ( u k b ) − ψ ( u k a ) ] , Δ W s , i h = 1 ∣ Ω i ∣ I i [ W s h , b − W s h , a ] , \Delta W_{s,i}^{p}=\frac{1}{|\Omega_i|}\sum_{k\in\mathcal{P}_{s,i}}m_k\left[\psi(u_k^b)-\psi(u_k^a)\right],\qquad\Delta W_{s,i}^{h}=\frac{1}{|\Omega_i|}\mathcal{I}_i\left[W_s^{h,b}-W_s^{h,a}\right], Δ W s , i p = ∣ Ω i ∣ 1 k ∈ P s , i ∑ m k [ ψ ( u k b ) − ψ ( u k a ) ] , Δ W s , i h = ∣ Ω i ∣ 1 I i [ W s h , b − W s h , a ] , W s , i b = W s , i a + Δ W s , i p + Δ W s , i h , W s , i h , b = W s , i h , a + Δ W s , i h . (23) W_{s,i}^b=W_{s,i}^a+\Delta W_{s,i}^p+\Delta W_{s,i}^h,\qquad W_{s,i}^{h,b}=W_{s,i}^{h,a}+\Delta W_{s,i}^h. \tag*{(23)} W s , i b = W s , i a + Δ W s , i p + Δ W s , i h , W s , i h , b = W s , i h , a + Δ W s , i h . (23) As the moments (22) are polynomials of degree at most three in each coordinate, I i \mathcal{I}_i I i W s , i h , b W_{s,i}^{h,b} W s , i h , b h s b h_s^b h s b
The total energy of a state with the cell averages W s , i = ( ρ s , i , P s , i , E s , i ) T W_{s,i}=(\rho_{s,i},P_{s,i},\mathcal{E}_{s,i})^{\mathsf T} W s , i = ( ρ s , i , P s , i , E s , i ) T e e e b b b
E t o t = ∑ s , i ∣ Ω i ∣ E s , i + ϵ 0 2 e T M E e + 1 2 μ 0 b T M B b . \mathcal{E}_{\mathrm{tot}}=\sum_{s,i}|\Omega_i|\mathcal{E}_{s,i}+\frac{\epsilon_0}{2}e^{\mathsf T}M_Ee+\frac{1}{2\mu_0}b^{\mathsf T}M_Bb. E tot = s , i ∑ ∣ Ω i ∣ E s , i + 2 ϵ 0 e T M E e + 2 μ 0 1 b T M B b . Proposition 1. Let the rotation field b r o t b^{\mathrm{rot}} b rot
(i) Equation (18) has exactly one solution e ˉ \bar e e ˉ Δ t > 0 \Delta t>0 Δ t > 0
(ii) The particle and wave updates (14) and (20) –(23) change the energy of the species by ∑ s , i ∣ Ω i ∣ ( E s , i b − E s , i a ) = Δ t e ˉ T j ˉ \sum_{s,i}|\Omega_i|(\mathcal{E}_{s,i}^b-\mathcal{E}_{s,i}^a)=\Delta t\,\bar e^{\mathsf T}\bar j ∑ s , i ∣ Ω i ∣ ( E s , i b − E s , i a ) = Δ t e ˉ T j ˉ j ˉ \bar j j ˉ (16) .
(iii) With periodic or perfectly conducting field boundaries, ( W a , e n , b n ) (W^a,e^n,b^n) ( W a , e n , b n ) ( W b , e b , b b ) (W^b,e^b,b^b) ( W b , e b , b b )
Proof. (i) As L s L_s L s α s − 1 + α s − T = 2 I \alpha_s^{-1}+\alpha_s^{-\mathsf T}=2I α s − 1 + α s − T = 2 I z T α s z = ∣ α s z ∣ 2 z^{\mathsf T}\alpha_sz=|\alpha_sz|^2 z T α s z = ∣ α s z ∣ 2 (15) , (16) and ρ s h , a ≥ 0 \rho_s^{h,a}\ge0 ρ s h , a ≥ 0
w T X w = ∑ s κ s 2 ( ∑ k ∈ P s m k ∣ α k N k w ∣ 2 + ∑ i I i [ ρ s h , a ∣ α s N E w ∣ 2 ] ) ≥ 0 w^{\mathsf T}Xw=\sum_s\kappa_s^2\left(\sum_{k\in\mathcal{P}_s}m_k|\alpha_kN_kw|^2+\sum_i\mathcal{I}_i\left[\rho_s^{h,a}|\alpha_sN_Ew|^2\right]\right)\ge0 w T X w = s ∑ κ s 2 ( k ∈ P s ∑ m k ∣ α k N k w ∣ 2 + i ∑ I i [ ρ s h , a ∣ α s N E w ∣ 2 ] ) ≥ 0 for every w w w M E M_E M E C T M B − 1 C C^{\mathsf T}M_B^{-1}C C T M B − 1 C (18) is at least ϵ 0 w T M E w > 0 \epsilon_0w^{\mathsf T}M_Ew>0 ϵ 0 w T M E w > 0 w ≠ 0 w\ne0 w = 0
(ii) For a particle, (13) , (14) and u ˉ k ⋅ ( u ˉ k × B r o t ) = 0 \bar u_k\cdot(\bar u_k\times B^{\mathrm{rot}})=0 u ˉ k ⋅ ( u ˉ k × B rot ) = 0
1 2 ( ∣ u k b ∣ 2 − ∣ u k a ∣ 2 ) = u ˉ k ⋅ ( u k b − u k a ) = Δ t κ s e ˉ T N k T α k ( u k a + δ κ s N k e ˉ ) , \frac{1}{2}\left(|u_k^b|^2-|u_k^a|^2\right)=\bar u_k\cdot(u_k^b-u_k^a)=\Delta t\,\kappa_s\bar e^{\mathsf T}N_k^{\mathsf T}\alpha_k\left(u_k^a+\delta\kappa_sN_k\bar e\right), 2 1 ( ∣ u k b ∣ 2 − ∣ u k a ∣ 2 ) = u ˉ k ⋅ ( u k b − u k a ) = Δ t κ s e ˉ T N k T α k ( u k a + δ κ s N k e ˉ ) , whose sum with the weights m k m_k m k Δ t e ˉ T ( j ^ p + δ X p e ˉ ) \Delta t\,\bar e^{\mathsf T}(\hat j_p+\delta X_p\bar e) Δ t e ˉ T ( j ^ p + δ X p e ˉ ) (15) . The wave is pushed by the same velocity map (20) as the particles, so the identity holds at every point x x x u u u h s a ( x , u ) d u h_s^a(x,u)\,\mathrm{d}u h s a ( x , u ) d u (21) , E s h , b − E s h , a = Δ t κ s e ˉ T N E T P ˉ s h \mathcal{E}_s^{h,b}-\mathcal{E}_s^{h,a}=\Delta t\,\kappa_s\bar e^{\mathsf T}N_E^{\mathsf T}\bar P_s^h E s h , b − E s h , a = Δ t κ s e ˉ T N E T P ˉ s h (23) sum by (16) to Δ t e ˉ T ( j ^ h + δ X h e ˉ ) \Delta t\,\bar e^{\mathsf T}(\hat j_h+\delta X_h\bar e) Δ t e ˉ T ( j ^ h + δ X h e ˉ )
(iii) The fields (19) satisfy (17) with b ˉ = 1 2 ( b b + b n ) = b n − δ M B − 1 C e ˉ \bar b=\frac{1}{2}(b^b+b^n)=b^n-\delta M_B^{-1}C\bar e b ˉ = 2 1 ( b b + b n ) = b n − δ M B − 1 C e ˉ M E M_E M E M B M_B M B M B ( b b − b n ) = − Δ t C e ˉ M_B(b^b-b^n)=-\Delta t\,C\bar e M B ( b b − b n ) = − Δ t C e ˉ ϵ 0 e ˉ T M E ( e b − e n ) + μ 0 − 1 b ˉ T M B ( b b − b n ) = Δ t μ 0 − 1 ( e ˉ T C T b ˉ − b ˉ T C e ˉ ) − Δ t e ˉ T j ˉ = − Δ t e ˉ T j ˉ \epsilon_0\bar e^{\mathsf T}M_E(e^b-e^n)+\mu_0^{-1}\bar b^{\mathsf T}M_B(b^b-b^n)=\Delta t\,\mu_0^{-1}(\bar e^{\mathsf T}C^{\mathsf T}\bar b-\bar b^{\mathsf T}C\bar e)-\Delta t\,\bar e^{\mathsf T}\bar j=-\Delta t\,\bar e^{\mathsf T}\bar j ϵ 0 e ˉ T M E ( e b − e n ) + μ 0 − 1 b ˉ T M B ( b b − b n ) = Δ t μ 0 − 1 ( e ˉ T C T b ˉ − b ˉ T C e ˉ ) − Δ t e ˉ T j ˉ = − Δ t e ˉ T j ˉ
Proposition 1 holds for any rotation field, so that the energy balance does not constrain its choice. The implicit midpoint rule evaluates the Lorentz force at the midpoint t n + δ t^{n+\delta} t n + δ t n + δ t^{n+\delta} t n + δ S 2 e n \mathbf{S}^{2}e^n S 2 e n ω p , i Δ t \omega_{p,i}\Delta t ω p , i Δ t β i \beta_i β i
b r o t = b n − δ Φ M B − 1 C S 2 e n , (24) \mathbf{b}^{\mathrm{rot}}=\mathbf{b}^n-\delta\mathbf{\Phi}\mathbf{M}_B^{-1}\mathbf{C}\mathbf{S}^{2}e^n, \tag*{(24)} b rot = b n − δ Φ M B − 1 C S 2 e n , (24) where S \mathbf{S} S e e e ( 1 4 , 1 2 , 1 4 ) \left(\frac{1}{4},\frac{1}{2},\frac{1}{4}\right) ( 4 1 , 2 1 , 4 1 ) e n e^n e n b r o t \mathbf{b}^{\mathrm{rot}} b rot Φ \mathbf{\Phi} Φ Ω i \Omega_i Ω i
Φ i = ( 1 + ( ω p , i Δ t ) 2 48 + β i 2 16 ) − 1 , β i = k ∗ V i ω p , i Δ t 2 , k ∗ = 2 h min , \Phi_i=\left(1+\frac{(\omega_{p,i}\Delta t)^2}{48}+\frac{\beta_i^2}{16}\right)^{-1},\qquad\beta_i=k_*V_i\omega_{p,i}\Delta t^2,\qquad k_*=\frac{2}{h_{\min}}, Φ i = ( 1 + 48 ( ω p , i Δ t ) 2 + 16 β i 2 ) − 1 , β i = k ∗ V i ω p , i Δ t 2 , k ∗ = h m i n 2 , ω p , i 2 = 1 ϵ 0 ∑ s κ s 2 ρ s , i a , ω p , i 2 V i 2 = 1 ϵ 0 ∑ s κ s 2 ∣ P s , i a ∣ 2 ρ s , i a . \omega_{p,i}^2=\frac{1}{\epsilon_0}\sum_s\kappa_s^2\rho^a_{s,i},\qquad\omega_{p,i}^2V_i^2=\frac{1}{\epsilon_0}\sum_s\kappa_s^2\frac{\lvert\mathbf{P}^a_{s,i}\rvert^2}{\rho^a_{s,i}}. ω p , i 2 = ϵ 0 1 s ∑ κ s 2 ρ s , i a , ω p , i 2 V i 2 = ϵ 0 1 s ∑ κ s 2 ρ s , i a ∣ P s , i a ∣ 2 . Here ω p , i \omega_{p,i} ω p , i V i V_i V i Ω i \Omega_i Ω i h min h_{\min} h m i n k ∗ k_* k ∗ Φ i \Phi_i Φ i Δ t \Delta t Δ t ω p , i Δ t \omega_{p,i}\Delta t ω p , i Δ t sinc 2 ( ω p , i Δ t / 4 ) \operatorname{sinc}^2(\omega_{p,i}\Delta t/4) sinc 2 ( ω p , i Δ t /4 ) ( ω p , i Δ t ) − 2 (\omega_{p,i}\Delta t)^{-2} ( ω p , i Δ t ) − 2 b r o t \mathbf{b}^{\mathrm{rot}} b rot b n \mathbf{b}^n b n β i = 4 \beta_i=4 β i = 4 Φ \mathbf{\Phi} Φ S \mathbf{S} S (24) keeps a positive definite quadratic invariant that bounds the perturbations of a drifting cold plasma (Appendix § 7 ).
Gauss correction
The Gauss correction G \mathcal{G} G (28) . It leaves the positions and velocities of the particles and the magnetic field unchanged.
With the moments W s h , c ( x ) = ( ρ s h , c , P s h , c , E s h , c ) T \mathbf{W}^{h,c}_s(\mathbf{x})=(\rho^{h,c}_s,\mathbf{P}^{h,c}_s,\mathcal{E}^{h,c}_s)^{\mathsf{T}} W s h , c ( x ) = ( ρ s h , c , P s h , c , E s h , c ) T h s c h^c_s h s c
Q = ∑ s κ s ( ∑ k ∈ P s m k Λ ( x k ) T + ∫ Ω Λ T ρ s h , c d x ) , D = ∫ Ω Λ T ( ∇ ⋅ N E ) d x , (25) \mathbf{Q}=\sum_s\kappa_s\left(\sum_{k\in\mathcal{P}_s}m_k\mathbf{\Lambda}(\mathbf{x}_k)^{\mathsf{T}}+\int_\Omega\mathbf{\Lambda}^{\mathsf{T}}\rho^{h,c}_s\,d\mathbf{x}\right),\qquad\mathbf{D}=\int_\Omega\mathbf{\Lambda}^{\mathsf{T}}(\nabla\cdot\mathbf{N}_E)\,d\mathbf{x}, \tag*{(25)} Q = s ∑ κ s ( k ∈ P s ∑ m k Λ ( x k ) T + ∫ Ω Λ T ρ s h , c d x ) , D = ∫ Ω Λ T ( ∇ ⋅ N E ) d x , (25) and the Gauss defect is r = ϵ 0 D e c − Q \mathbf{r}=\epsilon_0\mathbf{D}e^c-\mathbf{Q} r = ϵ 0 D e c − Q r ℓ r_\ell r ℓ ϵ 0 ∇ ⋅ E h = ρ c \epsilon_0\nabla\cdot\mathbf{E}_h=\rho_c ϵ 0 ∇ ⋅ E h = ρ c λ ℓ \lambda_\ell λ ℓ ρ c \rho_c ρ c Δ e \Delta e Δ e r \mathbf{r} r ϵ 0 D Δ e \epsilon_0\mathbf{D}\Delta e ϵ 0 D Δ e D \mathbf{D} D G \mathcal{G} G r ∥ \mathbf{r}_{\parallel} r ∥ D \mathbf{D} D r ⊥ = r − r ∥ \mathbf{r}_{\perp}=\mathbf{r}-\mathbf{r}_{\parallel} r ⊥ = r − r ∥ D T r ⊥ = 0 \mathbf{D}^{\mathsf{T}}\mathbf{r}_{\perp}=0 D T r ⊥ = 0
For a charged species s s s Ω i \Omega_i Ω i W s , i c ≠ 0 W^c_{s,i}\ne0 W s , i c = 0 P s , i \mathcal{P}_{s,i} P s , i s s s Ω i \Omega_i Ω i 1 + ζ s 1+\zeta_s 1 + ζ s
ζ s ( x ) = ∑ ℓ ζ s , i , ℓ λ ℓ ( x ) ( x ∈ Ω i ) , ζ s , i , ℓ ≥ − 1 , (26) \zeta_s(\mathbf{x})=\sum_\ell\zeta_{s,i,\ell}\lambda_\ell(\mathbf{x})\quad(\mathbf{x}\in\Omega_i),\qquad\zeta_{s,i,\ell}\ge-1, \tag*{(26)} ζ s ( x ) = ℓ ∑ ζ s , i , ℓ λ ℓ ( x ) ( x ∈ Ω i ) , ζ s , i , ℓ ≥ − 1 , (26) where the sum runs over the nodes ℓ \ell ℓ Ω i \Omega_i Ω i λ ℓ ≥ 0 \lambda_\ell\ge0 λ ℓ ≥ 0 ∑ ℓ λ ℓ = 1 \sum_\ell\lambda_\ell=1 ∑ ℓ λ ℓ = 1 1 + ζ s ≥ 0 1+\zeta_s\ge0 1 + ζ s ≥ 0
m k n + 1 = ( 1 + ζ s ( x k ) ) m k , h s n + 1 ( x , u ) = ( 1 + ζ s ( x ) ) h s c ( x , u ) , m_k^{n+1} = \left(1+\zeta_s(x_k)\right)m_k,\qquad h_s^{n+1}(x,u)=\left(1+\zeta_s(x)\right)h_s^c(x,u), m k n + 1 = ( 1 + ζ s ( x k ) ) m k , h s n + 1 ( x , u ) = ( 1 + ζ s ( x ) ) h s c ( x , u ) , W s , i n + 1 = W s , i c + 1 ∣ Ω i ∣ ( ∑ k ∈ P s , i m k ζ s ( x k ) ψ ( u k ) + ∫ Ω i ζ s W s h , c d x ) , (27) W_{s,i}^{n+1}=W_{s,i}^c+\frac{1}{|\Omega_i|}\left(\sum_{k\in P_{s,i}}m_k\zeta_s(x_k)\psi(u_k)+\int_{\Omega_i}\zeta_sW_s^{h,c}\,dx\right), \tag*{(27)} W s , i n + 1 = W s , i c + ∣ Ω i ∣ 1 k ∈ P s , i ∑ m k ζ s ( x k ) ψ ( u k ) + ∫ Ω i ζ s W s h , c d x , (27) W s , i h , n + 1 = W s , i h , c + 1 ∣ Ω i ∣ ∫ Ω i ζ s W s h , c d x , W_{s,i}^{h,n+1}=W_{s,i}^{h,c}+\frac{1}{|\Omega_i|}\int_{\Omega_i}\zeta_sW_s^{h,c}\,dx, W s , i h , n + 1 = W s , i h , c + ∣ Ω i ∣ 1 ∫ Ω i ζ s W s h , c d x , with the electric field e n + 1 e^{n+1} e n + 1 Q n + 1 Q^{n+1} Q n + 1 m k n + 1 m_k^{n+1} m k n + 1 ( 1 + ζ s ) ρ s h (1+\zeta_s)\rho_s^h ( 1 + ζ s ) ρ s h m k m_k m k ρ s h \rho_s^h ρ s h Q n + 1 Q^{n+1} Q n + 1 W s , i n + 1 W_{s,i}^{n+1} W s , i n + 1 W s , i h , n + 1 W_{s,i}^{h,n+1} W s , i h , n + 1
The field e n + 1 e^{n+1} e n + 1 ζ s , i , ℓ \zeta_{s,i,\ell} ζ s , i , ℓ
J = ϵ 0 2 ( e n + 1 − e c ) T M E ( e n + 1 − e c ) + 1 2 ∑ s , i k B T s , i m s ( ∑ k ∈ P s , i m k ζ s ( x k ) 2 + ∫ Ω i ρ s h , c ζ s 2 d x ) (28) \mathcal{J}=\frac{\epsilon_0}{2}(e^{n+1}-e^c)^{\mathsf{T}}M_E(e^{n+1}-e^c)+\frac{1}{2}\sum_{s,i}\frac{k_BT_{s,i}}{m_s}\left(\sum_{k\in P_{s,i}}m_k\zeta_s(x_k)^2+\int_{\Omega_i}\rho_s^{h,c}\zeta_s^2\,dx\right) \tag*{(28)} J = 2 ϵ 0 ( e n + 1 − e c ) T M E ( e n + 1 − e c ) + 2 1 s , i ∑ m s k B T s , i k ∈ P s , i ∑ m k ζ s ( x k ) 2 + ∫ Ω i ρ s h , c ζ s 2 d x (28) subject to the bounds in (26) and to
ϵ 0 D e n + 1 − Q n + 1 = r ⊥ , ∑ i ∣ Ω i ∣ ( W s , i n + 1 − W s , i c ) = 0 ( κ s ≠ 0 ) , (29) \epsilon_0De^{n+1}-Q^{n+1}=r_\perp,\qquad\sum_i|\Omega_i|\left(W_{s,i}^{n+1}-W_{s,i}^c\right)=0\quad(\kappa_s\ne0), \tag*{(29)} ϵ 0 D e n + 1 − Q n + 1 = r ⊥ , i ∑ ∣ Ω i ∣ ( W s , i n + 1 − W s , i c ) = 0 ( κ s = 0 ) , (29) so that the Gauss defect after the correction is r ⊥ r_\perp r ⊥ J \mathcal{J} J ( k B T s , i / m s ) ∫ ρ ln ρ d x (k_BT_{s,i}/m_s)\int\rho\ln\rho\,dx ( k B T s , i / m s ) ∫ ρ ln ρ d x ρ \rho ρ 1 + ζ s 1+\zeta_s 1 + ζ s
For a uniform plasma on a periodic domain, let r = ϵ 0 ∇ ⋅ E − ρ c r=\epsilon_0\nabla\cdot E-\rho_c r = ϵ 0 ∇ ⋅ E − ρ c Δ E = − ∇ ϕ \Delta E=-\nabla\phi Δ E = − ∇ ϕ ζ s = − κ s m s ϕ / ( k B T s ) \zeta_s=-\kappa_sm_s\phi/(k_BT_s) ζ s = − κ s m s ϕ / ( k B T s ) − ∇ 2 ϕ + λ D − 2 ϕ = − r / ϵ 0 -\nabla^2\phi+\lambda_D^{-2}\phi=-r/\epsilon_0 − ∇ 2 ϕ + λ D − 2 ϕ = − r / ϵ 0 λ D \lambda_D λ D ϕ \phi ϕ − r -r − r λ D − 2 ϕ \lambda_D^{-2}\phi λ D − 2 ϕ ϵ 0 → 0 \epsilon_0\to0 ϵ 0 → 0 1 / ϵ 0 1/\epsilon_0 1/ ϵ 0 J \mathcal{J} J ζ = 0 \zeta=0 ζ = 0 J \mathcal{J} J
The correction has the following properties, proved in Appendix § 8 .
Proposition 2.
(i) The minimization of (28) subject to (26) and (29) has a solution, and all solutions give the same corrected state (27).
(ii) In this state, the Gauss defect ϵ 0 D e n + 1 − Q n + 1 \epsilon_0De^{n+1}-Q^{n+1} ϵ 0 D e n + 1 − Q n + 1 r ⊥ r_\perp r ⊥
(iii) The total energy and the electric energy ϵ 0 2 e T M E e \frac{\epsilon_0}{2}e^{\mathsf{T}}M_Ee 2 ϵ 0 e T M E e
The minimization is a convex quadratic problem whose only inequality constraints are the bounds in (26). It is solved by the primal–dual active-set method [29 ] . Each iteration eliminates the coefficients cell by cell and solves one linear system for e n + 1 e^{n+1} e n + 1
Algorithm
One time step of the method from t n t^n t n t n + 1 t^{n+1} t n + 1 Figure 2 ). Steps 1 and 3 are the UGKWP step of § 3.2 as modified in § 4.2 .
Step 1. Perform the first transport step A δ \mathcal{A}_{\delta} A δ h s n h_s^n h s n τ s , i \tau_{s,i} τ s , i W s , i n \mathbf{W}_{s,i}^n W s , i n W s , i a \mathbf{W}_{s,i}^a W s , i a W s , i h , a \mathbf{W}_{s,i}^{h,a} W s , i h , a
Step 2. Perform the field step V Δ t \mathcal{V}_{\Delta t} V Δ t h s a h_s^a h s a W s , i h , a \mathbf{W}_{s,i}^{h,a} W s , i h , a e ˉ \bar{\mathbf{e}} e ˉ e b \mathbf{e}^b e b b b \mathbf{b}^b b b h s b h_s^b h s b W s , i b \mathbf{W}_{s,i}^b W s , i b W s , i h , b \mathbf{W}_{s,i}^{h,b} W s , i h , b
Step 3. Perform the second transport step A δ \mathcal{A}_{\delta} A δ τ s , i \tau_{s,i} τ s , i W s , i b \mathbf{W}_{s,i}^b W s , i b h s b h_s^b h s b h s n h_s^n h s n
Step 4. Apply the Gauss correction G \mathcal{G} G h s c h_s^c h s c W s , i h , c \mathbf{W}_{s,i}^{h,c} W s , i h , c § 4.4 , and set e n + 1 \mathbf{e}^{n+1} e n + 1 m k n + 1 m_k^{n+1} m k n + 1 h s n + 1 = ( 1 + ζ s ) h s c h_s^{n+1}=(1+\zeta_s)h_s^c h s n + 1 = ( 1 + ζ s ) h s c W s , i n + 1 \mathbf{W}_{s,i}^{n+1} W s , i n + 1 W s , i h , n + 1 \mathbf{W}_{s,i}^{h,n+1} W s , i h , n + 1 b n + 1 = b b \mathbf{b}^{n+1}=\mathbf{b}^b b n + 1 = b b
With periodic boundaries, the transport steps change neither e \mathbf{e} e b \mathbf{b} b
E t o t n + 1 − E t o t n = ϵ 0 2 [ ( e n + 1 ) T M E e n + 1 − ( e c ) T M E e c ] , \mathcal{E}_{\mathrm{tot}}^{n+1}-\mathcal{E}_{\mathrm{tot}}^n=\frac{\epsilon_0}{2}\left[(\mathbf{e}^{n+1})^\mathsf{T}M_E\mathbf{e}^{n+1}-(\mathbf{e}^c)^\mathsf{T}M_E\mathbf{e}^c\right], E tot n + 1 − E tot n = 2 ϵ 0 [ ( e n + 1 ) T M E e n + 1 − ( e c ) T M E e c ] , which is the change of the electric energy in the Gauss correction.
Numerical tests
All quantities are nondimensional in the sense of § 2 . The electron plasma and cyclotron frequencies of the initial state are ω p e = ( κ e 2 ρ e / ϵ 0 ) 1 / 2 \omega_{pe}=(\kappa_e^2\rho_e/\epsilon_0)^{1/2} ω p e = ( κ e 2 ρ e / ϵ 0 ) 1/2 Ω c e = ∣ κ e ∣ ∣ B ∣ \Omega_{ce}=|\kappa_e||\mathbf{B}| Ω ce = ∣ κ e ∣∣ B ∣ y y y Δ y \Delta y Δ y
The time step is determined by the CFL condition on the flow and sound speeds of the species,
Δ t = C F L h min max s , i ( ∣ U s , i ∣ + c s , i ) , c s , i = 5 k B T s , i 3 m s , (30) \Delta t=\mathrm{CFL}\,\frac{h_{\min}}{\max_{s,i}(|\mathbf{U}_{s,i}|+c_{s,i})},\qquad c_{s,i}=\sqrt{\frac{5k_{\mathrm{B}}T_{s,i}}{3m_s}}, \tag*{(30)} Δ t = CFL max s , i ( ∣ U s , i ∣ + c s , i ) h m i n , c s , i = 3 m s 5 k B T s , i , (30) where c s , i c_{s,i} c s , i ∣ U s , i ∣ |\mathbf{U}_{s,i}| ∣ U s , i ∣
In one dimension, the RME (18) and the linear systems of the Gauss correction are solved by LU factorization. In two dimensions, the RME is solved by GMRES [30 ] with a restricted additive Schwarz preconditioner [31 ] , and the linear systems of the Gauss correction by MINRES [32 ] with a Chebyshev polynomial preconditioner. Both are solved to the relative residual 10 − 8 10^{-8} 1 0 − 8 § 5.3 , and in § 5.4 the RME to 10 − 10 10^{-10} 1 0 − 10 10 − 6 10^{-6} 1 0 − 6 6 × 10 − 12 6\times10^{-12} 6 × 1 0 − 12
The comparison method UGKWP-PHM is the UGKWP method of [22 ] with the same collision model. It advances the perfectly hyperbolic Maxwell (PHM) equations [23 ] explicitly by a second-order finite-volume scheme with the monotonized central (MC) limiter. The damping rate of the cleaning fields is σ ϕ = 0 \sigma_{\phi}=0 σ ϕ = 0 σ ϕ = 3 \sigma_{\phi}=3 σ ϕ = 3 (30) with max ( ∣ U s , i ∣ + c s , i , c ) \max(\lvert U_{s,i}\rvert+c_{s,i},c) max (∣ U s , i ∣ + c s , i , c ) ∣ U s , i ∣ + c s , i \lvert U_{s,i}\rvert+c_{s,i} ∣ U s , i ∣ + c s , i = 0.5 =0.5 = 0.5 0.4 0.4 0.4 = 0.5 =0.5 = 0.5
The CPU times reported in § 5.3 and § 5.4 are wall-clock times multiplied by the number of cores. They are measured on 64 cores of a workstation with two 64-core AMD EPYC 7742 processors.
Plasma waves
Plasma waves test the longitudinal and the magnetized transverse response of the field step and its energy balance. An electron plasma wave covers the transition from collisionless Landau damping to collision-dominated oscillations, and a parallel whistler wave in a warm plasma tests the frequency shift from the cold-plasma value and the electron cyclotron damping. With a collision frequency that varies over four orders of magnitude in space, the electron plasma wave also tests the method at the transitions between weakly collisional and collision-dominated regions.
The electron plasma wave is computed with m e = 1 m_e=1 m e = 1 κ e = − 1 \kappa_e=-1 κ e = − 1 ϵ 0 = 1 \epsilon_0=1 ϵ 0 = 1 c = 20 c=20 c = 20 [ 0 , 2 π ] [0,2\pi] [ 0 , 2 π ] T e = 0.25 T_e=0.25 T e = 0.25
n e ( y , 0 ) = 1 − a k cos k y , k = 1 , a = 0.05 , n_e(y,0)=1-ak\cos ky,\qquad k=1,\qquad a=0.05, n e ( y , 0 ) = 1 − ak cos k y , k = 1 , a = 0.05 , so that ω p e = 1 \omega_{pe}=1 ω p e = 1 k λ D e = 0.5 k\lambda_{De}=0.5 k λ D e = 0.5 λ D e = ( T e / m e ) 1 / 2 / ω p e \lambda_{De}=(T_e/m_e)^{1/2}/\omega_{pe} λ D e = ( T e / m e ) 1/2 / ω p e m i = 1836 m_i=1836 m i = 1836 κ i = 1 / 1836 \kappa_i=1/1836 κ i = 1/1836 T i = 0.01 T_i=0.01 T i = 0.01 ν e = χ e e n e \nu_e=\chi_{ee}n_e ν e = χ ee n e χ e e = 0 \chi_{ee}=0 χ ee = 0 1 1 1 10 10 10 ν e = 0 \nu_e=0 ν e = 0 ω p e \omega_{pe} ω p e 10 ω p e 10\omega_{pe} 10 ω p e ν e = 0 \nu_e=0 ν e = 0 y 0 + a sin k y 0 y_0+a\sin ky_0 y 0 + a sin k y 0 y 0 y_0 y 0 E y E_y E y B = 0 B=0 B = 0 4 × 10 4 4\times10^4 4 × 1 0 4 ω \omega ω γ \gamma γ A ∞ + A 1 e γ t cos ( ω t + ϕ ) A_\infty+A_1e^{\gamma t}\cos(\omega t+\phi) A ∞ + A 1 e γ t cos ( ω t + ϕ ) A ( t ) A(t) A ( t ) sin k y \sin ky sin k y E y E_y E y 2 ≤ t ≤ 12 2\le t\le12 2 ≤ t ≤ 12 A ∞ A_\infty A ∞ ν e \nu_e ν e [ − 4.5 , 4.5 ] [-4.5,4.5] [ − 4.5 , 4.5 ]
As ν e \nu_e ν e ( 1 + 5 3 k 2 λ D e 2 ) 1 / 2 ω p e = 1.190 (1+\frac{5}{3}k^2\lambda_{De}^2)^{1/2}\omega_{pe}=1.190 ( 1 + 3 5 k 2 λ D e 2 ) 1/2 ω p e = 1.190 Figure 3 , Table 1 ). At ν e = 10 ω p e \nu_e=10\omega_{pe} ν e = 10 ω p e 1 − e − ν e Δ t / 2 ≈ 0.44 1-e^{-\nu_e\Delta t/2}\approx0.44 1 − e − ν e Δ t /2 ≈ 0.44 0.21 % 0.21\% 0.21% 0.009 ω p e 0.009\omega_{pe} 0.009 ω p e Table 1 ). At ν e = 10 ω p e \nu_e=10\omega_{pe} ν e = 10 ω p e A ( t ) A(t) A ( t ) Figure 3 (f)). For ν e = 0 \nu_e=0 ν e = 0 ω p e \omega_{pe} ω p e c k = 20 ck=20 c k = 20 Figure 3 (d), (Figure 3 e)). It is the oscillation of the cleaning field of the perfectly hyperbolic Maxwell equations, which travels at the speed of light and is undamped for σ ϕ = 0 \sigma_\phi= 0 σ ϕ = 0
Figure 3. Electron plasma wave with (a), (d) ν e = 0 \nu_e = 0 ν e = 0 ν e = ω p e \nu_e = \omega_{pe} ν e = ω p e ν e = 10 ω p e \nu_e = 10\omega_{pe} ν e = 10 ω p e ∣ A ( t ) ∣ / A ( 0 ) |A(t)|/A(0) ∣ A ( t ) ∣/ A ( 0 ) e γ L t e^{\gamma_L t} e γ L t A ( t ) A(t) A ( t ) sin k y \sin ky sin k y E y E_y E y γ L = − 0.1533 \gamma_L = -0.1533 γ L = − 0.1533 A ( t ) A(t) A ( t ) A ( 0 ) A(0) A ( 0 )
Table 1. Plasma waves: frequency ω \omega ω γ \gamma γ
The parallel whistler wave propagates along the uniform magnetic field B 0 = 0.25 y ^ \mathbf{B}_0 = 0.25\hat{\mathbf{y}} B 0 = 0.25 y ^ m i / m e = 100 m_i/m_e = 100 m i / m e = 100 T e = T i = 1.25 T_e = T_i = 1.25 T e = T i = 1.25 ω p e / Ω c e = 4 \omega_{pe}/\Omega_{ce} = 4 ω p e / Ω ce = 4 2 π / k 2\pi/k 2 π / k k d e = 1.25 k d_e = 1.25 k d e = 1.25 d e = c / ω p e d_e = c/\omega_{pe} d e = c / ω p e ∣ B ⊥ ∣ = 0.01 B 0 |B_\perp| = 0.01B_0 ∣ B ⊥ ∣ = 0.01 B 0 ω \omega ω γ \gamma γ A 1 e − i ( ω + i γ ) t A_1 e^{-i(\omega+i\gamma)t} A 1 e − i ( ω + iγ ) t B x + i B z B_x + iB_z B x + i B z k k k 10 ≤ t ≤ 80 10 \le t \le80 10 ≤ t ≤ 80 ω = 0.1519 \omega= 0.1519 ω = 0.1519 B 0 B_0 B 0
the electrons in cyclotron resonance with the wave absorb its energy. With a time step of about 11 light-crossing times of a cell, UGKWP-EC reproduces both effects (Table 1 ).
Without the Gauss correction, the computations of the three electron plasma waves and of the whistler wave conserve the total energy E t o t \mathcal{E}_{\mathrm{tot}} E tot § 4.3 to rounding error, with particles alone and with particles and wave (Proposition 1(iii), Figure 4 ). With the correction, the change of E t o t \mathcal{E}_{\mathrm{tot}} E tot G \mathcal{G} G 3 × 10 − 6 3 \times10^{-6} 3 × 1 0 − 6 E t o t \mathcal{E}_{\mathrm{tot}} E tot 10 − 4 10^{-4} 1 0 − 4
Figure 4. Relative change ∣ Δ E t o t ∣ / E t o t ( 0 ) |\Delta\mathcal{E}_{\mathrm{tot}}|/\mathcal{E}_{\mathrm{tot}}(0) ∣Δ E tot ∣/ E tot ( 0 ) Δ E t o t = E t o t ( t ) − E t o t ( 0 ) \Delta\mathcal{E}_{\mathrm{tot}}=\mathcal{E}_{\mathrm{tot}}(t)-\mathcal{E}_{\mathrm{tot}}(0) Δ E tot = E tot ( t ) − E tot ( 0 ) ν e = 0 \nu_e=0 ν e = 0 ν e = 10 ω p e \nu_e=10\omega_{pe} ν e = 10 ω p e G \mathcal{G} G
Finally, the electron plasma wave is computed on 64 and 128 cells with the electron collision frequency ν e = χ e e ( y ) n e \nu_e=\chi_{ee}(y)n_e ν e = χ ee ( y ) n e
χ e e ( y ) = 0.1 × 10 4 σ ( y ) , σ ( y ) = 1 2 + tanh ( 3 sin y ) 2 tanh 3 , (31) \chi_{ee}(y)=0.1 \times10^{4\sigma(y)}, \qquad\sigma(y)=\frac{1}{2}+\frac{\tanh(3\sin y)}{2\tanh3}, \tag*{(31)} χ ee ( y ) = 0.1 × 1 0 4 σ ( y ) , σ ( y ) = 2 1 + 2 tanh 3 tanh ( 3 sin y ) , (31) where χ e e \chi_{ee} χ ee y = 3 π / 2 y=3\pi/2 y = 3 π /2 y = π / 2 y=\pi/2 y = π /2 y = 0 y=0 y = 0 y = π y=\pi y = π Figure 5 (b)). The reference is DVM-Poisson with the same χ e e ( y ) \chi_{ee}(y) χ ee ( y )
Figure 5. Electron plasma wave with the electron collision coefficient (31): (a) E y E_y E y t = 5 t=5 t = 5 ρ e p / ρ e \rho^p_e/\rho_e ρ e p / ρ e t ≈ 2.5 t \approx2.5 t ≈ 2.5 χ e e ( y ) \chi_{ee}(y) χ ee ( y ) E y E_y E y ϖ e , i = 1 \varpi_{e,i}=1 ϖ e , i = 1
On 64 cells, ν e Δ t \nu_e\Delta t ν e Δ t Figure 5 (b)). The electric field E y E_y E y Figure 5 (a), (Figure 5 c)).
At the transitions, where the particle fraction varies in space, the analytic fluxes bring momentum into the cell averages but not into the particles, and this momentum is assigned to the wave through the residual momentum in (12). Without the residual momentum in the wave, that is, with ϖ e , i = 1 \varpi_{e,i}=1 ϖ e , i = 1 E y E_y E y Figure 5 (a), (Figure 5 c)).
Multiscale Brio–Wu problem
The Brio–Wu problem [33 ] for electrons and ions [34 ] tests whether the method follows kinetic solutions at finite Larmor radius and mean free path and approaches ideal MHD as both decrease. In this subsection, x x x m e = 0.01 m_e=0.01 m e = 0.01 m i = 1 m_i=1 m i = 1
( n s , T s , B ) = { ( 1 , 0.5 , ( 0.75 , 1 , 0 ) ) , x < 0.5 , ( 0.125 , 0.4 , ( 0.75 , − 1 , 0 ) ) , x > 0.5 , E = 0. (n_s,T_s,\mathbf{B}) =
\begin{cases}
(1,0.5,(0.75,1,0)), & x < 0.5, \\
(0.125,0.4,(0.75,-1,0)), & x > 0.5,
\end{cases}
\qquad\mathbf{E} = 0. ( n s , T s , B ) = { ( 1 , 0.5 , ( 0.75 , 1 , 0 )) , ( 0.125 , 0.4 , ( 0.75 , − 1 , 0 )) , x < 0.5 , x > 0.5 , E = 0. In the normalization of [8 ] , with the Knudsen number K n \mathrm{Kn} Kn r r r λ \lambda λ
κ i = 1 r , κ e = − 1 m e r , ϵ 0 = λ 2 , μ 0 − 1 = λ 2 c 2 , χ i i = 1 K n , χ e e = 1 K n m e . \kappa_i = \frac{1}{r}, \qquad\kappa_e = -\frac{1}{m_e r}, \qquad\epsilon_0 = \lambda^2, \qquad\mu_0^{-1} = \lambda^2 c^2, \qquad\chi_{ii} = \frac{1}{\mathrm{Kn}}, \qquad\chi_{ee} = \frac{1}{\mathrm{Kn}\sqrt{m_e}}. κ i = r 1 , κ e = − m e r 1 , ϵ 0 = λ 2 , μ 0 − 1 = λ 2 c 2 , χ ii = Kn 1 , χ ee = Kn m e 1 . with λ = 0.01 \lambda=0.01 λ = 0.01 c = 100 c=100 c = 100 ( r , K n ) = ( 0.1 , 1 ) (r,\mathrm{Kn})=(0.1,1) ( r , Kn ) = ( 0.1 , 1 ) ( 0.01 , 0.01 ) (0.01,0.01) ( 0.01 , 0.01 ) ( 10 − 3 , 10 − 3 ) (10^{-3},10^{-3}) ( 1 0 − 3 , 1 0 − 3 ) ω p e Δ t \omega_{pe}\Delta t ω p e Δ t Ω c e Δ t \Omega_{ce}\Delta t Ω ce Δ t c Δ t / Δ x c\Delta t/\Delta x c Δ t /Δ x [ 0 , 1 ] [0,1] [ 0 , 1 ] t = 0.12 t=0.12 t = 0.12
The first two settings are compared with kinetic solutions of iPIC3D [12 ] extended by a BGK operator (iPIC3D-BGK), and the third with ideal MHD. After every push, iPIC3D-BGK replaces each particle by a sample of the Maxwellian of its species in its cell with the probability 1 − e − ν s Δ t 1-e^{-\nu_s\Delta t} 1 − e − ν s Δ t ν s = χ s s n s \nu_s=\chi_{ss}n_s ν s = χ ss n s
iPIC3D-BGK computations use the same domain, a 10 times finer mesh, a 22 times smaller time step and 1024 particles per species and cell, and their profiles are averaged onto the 400 cells. The ideal-MHD solution is computed by the HLLD solver [35 ] on the same 400 cells, with γ = 5 / 3 \gamma= 5/3 γ = 5/3 m i n i + m e n e m_i n_i + m_e n_e m i n i + m e n e n i T i + n e T e n_i T_i + n_e T_e n i T i + n e T e
As r r r Figure 6 ). Like UGKWP-PHM, UGKWP-EC reproduces the profiles of iPIC3D-BGK for r = 0.1 r = 0.1 r = 0.1 Figure 6 (a), (Figure 6 b)). For r = 10 − 3 r = 10^{-3} r = 1 0 − 3 Figure 6 (c)). The pressure p = n i T i + n e T e p = n_i T_i + n_e T_e p = n i T i + n e T e Figure 6 (cFigure 5 ).
Figure 6. Brio–Wu problem at t = 0.12 t = 0.12 t = 0.12 r = 0.1 r = 0.1 r = 0.1 K n = 1 \mathrm{Kn} = 1 Kn = 1 r = 0.01 r = 0.01 r = 0.01 K n = 0.01 \mathrm{Kn} = 0.01 Kn = 0.01 r = 10 − 3 r = 10^{-3} r = 1 0 − 3 K n = 10 − 3 \mathrm{Kn} = 10^{-3} Kn = 1 0 − 3 u i u_i u i v i v_i v i x x x y y y B y B_y B y T i T_i T i p = n i T i + n e T e p = n_iT_i + n_eT_e p = n i T i + n e T e
In the MHD limit, the relative residual of Gauss’s law stays below 10 − 3 10^{-3} 1 0 − 3 10 − 2 10^{-2} 1 0 − 2
Orszag–Tang vortex
The Orszag–Tang vortex [36 ] tests the method in two dimensions, on a flow that forms shocks and thin current sheets, and shows how the number of particles and the computational cost change with the Knudsen number. It is computed on the periodic domain [ 0 , 2 π ] 2 [0, 2\pi]^2 [ 0 , 2 π ] 2 § 5.2 with m i = 1 m_i = 1 m i = 1 m e = 0.04 m_e = 0.04 m e = 0.04 λ = 0.05 \lambda= 0.05 λ = 0.05 c = 20 c = 20 c = 20 r = 0.2 r = 0.2 r = 0.2
ρ = γ 2 , p = γ , U = ( − sin y , sin x , 0 ) , B = ( − sin y , sin 2 x , 0 ) , E = − U × B , \rho= \gamma^2,\qquad p = \gamma,\qquad U = (-\sin y,\sin x,0),\qquad B = (-\sin y,\sin2x,0),\qquad E = -U \times B, ρ = γ 2 , p = γ , U = ( − sin y , sin x , 0 ) , B = ( − sin y , sin 2 x , 0 ) , E = − U × B , where γ = 5 / 3 \gamma= 5/3 γ = 5/3 ρ = ρ i + ρ e \rho= \rho_i + \rho_e ρ = ρ i + ρ e p = n i T i + n e T e p = n_i T_i + n_e T_e p = n i T i + n e T e U U U n i = n e n_i = n_e n i = n e T i = T e T_i = T_e T i = T e U U U
ρ i U i , z + ρ e U e , z = 0 , κ i ρ i U i , z + κ e ρ e U e , z = ( ∇ × B ) z μ 0 . \rho_i U_{i,z} + \rho_e U_{e,z} = 0,\qquad\kappa_i \rho_i U_{i,z} + \kappa_e \rho_e U_{e,z} = \frac{(\nabla\times B)_z}{\mu_0}. ρ i U i , z + ρ e U e , z = 0 , κ i ρ i U i , z + κ e ρ e U e , z = μ 0 ( ∇ × B ) z . The species have the collision coefficients χ i i = 1 / K n \chi_{ii} = 1/\mathrm{Kn} χ ii = 1/ Kn χ e e = 1 / ( K n m e ) \chi_{ee} = 1/(\mathrm{Kn}\sqrt{m_e}) χ ee = 1/ ( Kn m e ) K n = 10 − 3 \mathrm{Kn} = 10^{-3} Kn = 1 0 − 3 10 − 2 10^{-2} 1 0 − 2 10 − 1 10^{-1} 1 0 − 1 256 2 256^2 25 6 2 h h h t = 3 t = 3 t = 3 K n = 10 − 3 \mathrm{Kn} = 10^{-3} Kn = 1 0 − 3 10 − 1 10^{-1} 1 0 − 1 256 2 256^2 25 6 2 K n = 10 − 3 \mathrm{Kn} = 10^{-3} Kn = 1 0 − 3 h / 100 h/100 h /100 [37 ] with the finite-volume wave-propagation scheme of [38 ] and the locally implicit source update of [39 ] . For K n = 10 − 1 \mathrm{Kn} = 10^{-1} Kn = 1 0 − 1 § 5.2 ) with 4096 particles per species and cell and a fixed time step of about half the mean step of UGKWP-EC. The implicit scheme of iPIC3D-BGK solves the field equations for the electric field at t + θ Δ t t + \theta\Delta t t + θ Δ t θ = 1 / 2 \theta= 1/2 θ = 1/2
For both Knudsen numbers, UGKWP-EC forms the shocks, the central current sheet and the magnetic islands of the reference (Figure 7 and Figure 8 ). At K n = 10 − 1 \mathrm{Kn} = 10^{-1} Kn = 1 0 − 1 h h h K n = 10 − 3 \mathrm{Kn} = 10^{-3} Kn = 1 0 − 3
(a) Gkeyll five-moment
(b) Gkeyll five-moment
(c) UGKWP-EC
(d) UGKWP-EC
(e) UGKWP-PHM, CFL 0.4
(f) UGKWP-PHM, CFL 0.4
(g) UGKWP-PHM, CFL 0.5
(h) UGKWP-PHM, CFL 0.5
Figure 7. Orszag–Tang vortex, K n = 10 − 3 \mathrm{Kn}=10^{-3} Kn = 1 0 − 3 256 2 256^2 25 6 2 t = 3 t=3 t = 3 J z = ( ∇ × B ) z / μ 0 J_z=(\nabla\times\mathbf{B})_z/\mu_0 J z = ( ∇ × B ) z / μ 0 A z A_z A z J z J_z J z ∣ J z ∣ = 20 |J_z|=20 ∣ J z ∣ = 20
Figure 8. Orszag–Tang vortex, K n = 10 − 1 \mathrm{Kn}=10^{-1} Kn = 1 0 − 1 256 2 256^2 25 6 2 t = 3 t=3 t = 3 J z = ( ∇ × B ) z / μ 0 J_z=(\nabla\times B)_z/\mu_0 J z = ( ∇ × B ) z / μ 0 A z A_z A z J z J_z J z ∣ J z ∣ = 20 |J_z|=20 ∣ J z ∣ = 20
numbers, the pressure of UGKWP-EC along y = 0.625 π y = 0.625\pi y = 0.625 π 128 2 128^2 12 8 2 256 2 256^2 25 6 2 Figure 9 ).
Figure 9. Orszag–Tang vortex, t = 3 t = 3 t = 3 p p p y = 0.625 π y = 0.625\pi y = 0.625 π K n = 10 − 3 \mathrm{Kn} = 10^{-3} Kn = 1 0 − 3 K n = 10 − 1 \mathrm{Kn} = 10^{-1} Kn = 1 0 − 1
The total energy of UGKWP-EC changes only by the field-energy change of G \mathcal{G} G t = 3 t = 3 t = 3 Figure 7g , (Figure 7h ) and Figure 8 (g), (Figure 8 h)). At K n = 10 − 1 \mathrm{Kn} = 10^{-1} Kn = 1 0 − 1 t = 3 t = 3 t = 3 Figure 9 (b)). In the solution of this problem, the total thermal energy of each species does not fall below its initial value, because the entropy of each species does not decrease and the initial Maxwellian of uniform density and temperature has the largest entropy for its total mass and thermal energy. At t = 3 t = 3 t = 3
At K n = 10 − 3 \mathrm{Kn} = 10^{-3} Kn = 1 0 − 3 K n = 10 − 1 \mathrm{Kn} = 10^{-1} Kn = 1 0 − 1
The CPU time of UGKWP-EC is compared at the three Knudsen numbers with that of iPIC3D-BGK, which represents both species by particles at every Knudsen number. The iPIC3D-BGK computations use the same mesh, 256 particles per species and cell and θ = 1 \theta= 1 θ = 1 θ = 1 / 2 \theta= 1/2 θ = 1/2 ω p e Δ t ≈ 2.5 \omega_{\mathrm{pe}}\Delta t \approx2.5 ω pe Δ t ≈ 2.5 10 − 8 10^{-8} 1 0 − 8 Table 2 ). At CFL 0.4, UGKWP-PHM needs 4.0 to 5.7 times the CPU time of UGKWP-EC.
Table 2. Orszag–Tang vortex on 256 2 256^2 25 6 2
Magnetic reconnection
The collisionless Geospace Environmental Modeling (GEM) reconnection problem [40 ] tests the method on Hall reconnection, with time steps longer than the light-crossing time of a cell. It is computed on the domain [ − 12.8 , 12.8 ] × [ − 6.4 , 6.4 ] [-12.8, 12.8] \times[-6.4, 6.4] [ − 12.8 , 12.8 ] × [ − 6.4 , 6.4 ] m i / m e = 25 m_i/m_e = 25 m i / m e = 25 c / v A = 20 c/v_A = 20 c / v A = 20 v A v_A v A d i d_i d i Ω c i − 1 \Omega_{ci}^{-1} Ω c i − 1 B 0 B_0 B 0 n 0 n_0 n 0 v A = 1 v_A = 1 v A = 1 κ i = 1 \kappa_i = 1 κ i = 1 κ e = − 25 \kappa_e = -25 κ e = − 25 μ 0 = 1 \mu_0 = 1 μ 0 = 1 ϵ 0 = ( v A / c ) 2 \epsilon_0 = (v_A/c)^2 ϵ 0 = ( v A / c ) 2 T i = 5 / 12 T_i = 5/12 T i = 5/12 T e = 1 / 12 T_e = 1/12 T e = 1/12
n s s h = sech 2 ( 2 y ) , U s s h = U s z ^ , n s b g = 0.2 , U s b g = 0 , n_s^{\mathrm{sh}} = \operatorname{sech}^2(2y), \qquad U_s^{\mathrm{sh}} = U_s\hat{z}, \qquad n_s^{\mathrm{bg}} = 0.2, \qquad U_s^{\mathrm{bg}} = 0, n s sh = sech 2 ( 2 y ) , U s sh = U s z ^ , n s bg = 0.2 , U s bg = 0 , where U i = − 5 / 3 U_i = -5/3 U i = − 5/3 U e = 1 / 3 U_e = 1/3 U e = 1/3 E = 0 E = 0 E = 0 B = ∇ × ( A z z ^ ) B = \nabla\times(A_z\hat{z}) B = ∇ × ( A z z ^ )
A z = 0.5 ln cosh ( 2 y ) + 0.1 cos π x 12.8 cos π y 12.8 + 0.08 e − ( x 2 + y 2 ) / 4 , A_z = 0.5\ln\cosh(2y) + 0.1\cos\frac{\pi x}{12.8}\cos\frac{\pi y}{12.8} + 0.08e^{-(x^2+y^2)/4}, A z = 0.5 ln cosh ( 2 y ) + 0.1 cos 12.8 π x cos 12.8 π y + 0.08 e − ( x 2 + y 2 ) /4 , the sum of the Harris field, the GEM perturbation and a Gaussian perturbation at the X point. The particles carry the current of the Harris field only. The domain is periodic in x x x y = ± 6.4 y = \pm6.4 y = ± 6.4 256 × 128 256 \times128 256 × 128 h h h c / v A c/v_A c / v A 256 × 128 256 \times128 256 × 128 [8 ] on 256 × 128 256 \times128 256 × 128 40 3 40^3 4 0 3 [19 ] in velocity space. The reconnected flux Ψ = A z ( X ) − A z ( O ) \Psi= A_z(X) - A_z(O) Ψ = A z ( X ) − A z ( O ) A z A_z A z
At t = 12 t = 12 t = 12 Figure 10 ). The out-of-plane magnetic field of UGKWP-EC forms the Hall quadrupole of collisionless reconnection, with extremes along y = 0.5 y = 0.5 y = 0.5 Figure 11 (a)–(Figure 11 d), (Figure 11 f)). Through the X point, the current layer of UGKWP-EC is about as thin as in the references (Figure 11 (e)). Reconnection starts later in UGKWP-PHM at CFL 0.4, and more so at c / v A = 20 c/v_A = 20 c / v A = 20 Table 3 ), because the dissipation of its explicit Maxwell solver increases with c c c t = 12 t = 12 t = 12 Figure 11 ).
Figure 10. GEM problem, c / v A = 20 c/v_A = 20 c / v A = 20 256 × 128 256 \times128 256 × 128 J z = ( ∇ × B ) z J_z = (\nabla\times\mathbf{B})_z J z = ( ∇ × B ) z A z A_z A z ∣ y ∣ ≤ 4 |y| \le4 ∣ y ∣ ≤ 4 t = 12 t = 12 t = 12 24 24 24 ∣ J z ∣ = 2 |J_z| = 2 ∣ J z ∣ = 2
Table 3. GEM problem on 256 × 128 256 \times128 256 × 128 Ψ = 1 \Psi= 1 Ψ = 1 d Ψ / d t d\Psi/dt d Ψ/ d t Ψ ( 40 ) \Psi(40) Ψ ( 40 ) t = 40 t = 40 t = 40
Figure 11. GEM problem, c / v A = 20 c/v_A = 20 c / v A = 20 × \times × t = 12 t = 12 t = 12 B z B_z B z Ψ \Psi Ψ J z J_z J z y y y ∣ x ∣ ≤ 0.5 |x| \le0.5 ∣ x ∣ ≤ 0.5 B z B_z B z y = 0.5 y = 0.5 y = 0.5
For both values of c / v A c/v_A c / v A Figure 12 , Table 3 ). UGKWP-EC takes about 2200 steps for both values of c / v A c/v_A c / v A Table 4 ). UGKWP-PHM loses about 5% of its total energy, and at CFL 0.5 it reconnects still later and loses more than 19%, whereas the total energy of UGKWP-EC changes by less than 0.3% (Table 3 ).
Figure 12. GEM problem: reconnected flux Ψ ( t ) \Psi(t) Ψ ( t ) c / v A = 20 c/v_A = 20 c / v A = 20 c / v A = 10 c/v_A = 10 c / v A = 10 × \times ×
Table 4. GEM problem on 256 × 128 256 \times128 256 × 128
With G G G G G G Figure 13 (a)), and the total energy changes only by the field-energy change of G G G G G G Proposition 1 (iii), Figure 13 (b)), but the current layer through the X point breaks by t = 20 t = 20 t = 20 G \mathcal{G} G Figure 13 (c), (Figure 13 d)).
Figure 13. GEM problem, c / v A = 20 c/v_A = 20 c / v A = 20 256 × 128 256 \times128 256 × 128 κ i n 0 \kappa_i n_0 κ i n 0 G \mathcal{G} G Δ E G \Delta\mathcal{E}_{\mathcal{G}} Δ E G G \mathcal{G} G J z J_z J z t = 20 t = 20 t = 20 G \mathcal{G} G Figure 10 .
Conclusions
For the BGK–Vlasov–Maxwell system, this paper has presented UGKWP-EC, a semi-implicit UGKWP method in which the particles and the hydrodynamic wave share one electromagnetic response. The current response of the reformulated Maxwell equation is built from the implicit-midpoint velocity map that then pushes the particles and maps the wave. This equation is a single linear system for the electric field, with a unique solution for any time step. As Maxwell’s equations are advanced by the implicit midpoint rule with the same current, the field step conserves the discrete total energy exactly. The Gauss correction adjusts the electric field and the joint charge of particles and wave while conserving the mass, momentum and energy of every charged species. The transport steps sample and transport the field-pushed wave and the corrected wave without rebuilding them from their cell averages.
The field step remains stable without resolving the cyclotron period, the plasma period or the light-crossing time, so that the time step can be taken from a CFL condition on the flow and sound speeds of the species. The method reproduces the whistler wave of a warm plasma and the electron plasma wave from Landau damping to collision-dominated oscillations. For an electron plasma wave whose collision frequency varies over four orders of magnitude in space, it approaches the kinetic reference with mesh refinement. For the drifting plasma of Appendix A, its rotation field keeps the electrons as close to equilibrium as in the plasma without drift, whereas the rotation with the magnetic field of t n t^n t n ω p e Δ t ≈ 2.5 \omega_{pe}\Delta t \approx2.5 ω p e Δ t ≈ 2.5 10 − 3 10^{-3} 1 0 − 3 10 − 3 10^{-3} 1 0 − 3 c c c K n = 10 − 1 \mathrm{Kn}=10^{-1} Kn = 1 0 − 1 10 − 3 10^{-3} 1 0 − 3 40 3 40^3 4 0 3
Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
The current research is supported by the National Key R&D Program of China (Grant No. 2022YFA1004500), the National Natural Science Foundation of China (92371107), and the Hong Kong Research Grants Council (16208324).
Rotation field in a drifting plasma
The field step conserves the total energy for any rotation field (Proposition 1 ). For a drifting plasma, this does not bound the perturbations, because the drift can supply them with energy. This appendix first compares computations of a drifting plasma with the rotation field (24) and with the magnetic field B h n = N B n B_h^n=N B^n B h n = N B n t n t^n t n [13 ] uses. It then shows that, for a drifting cold plasma, the time step with (24), without Φ \boldsymbol{\Phi} Φ S \boldsymbol{S} S β i < 4 \beta_i < 4 β i < 4 § 4.3 . With the field of t n t^n t n
The plasma is the exact uniform equilibrium of collisionless electrons and ions with the density n 0 n_0 n 0 T T T
U e = U i = V = ( V , 0 , 0 ) , B = ( 0 , 0 , 1 ) , E = − V × B , U_e = U_i = V = (V,0,0), \qquad B = (0,0,1), \qquad E = -V \times B, U e = U i = V = ( V , 0 , 0 ) , B = ( 0 , 0 , 1 ) , E = − V × B , Figure A.1. Drifting plasma, T = 0.1 T = 0.1 T = 0.1 V = 0.5 V = 0.5 V = 0.5 3 2 T L \frac{3}{2}TL 2 3 T L G \mathcal{G} G Δ t = 0.025 \Delta t = 0.025 Δ t = 0.025 B h r o t = B h n \mathbf{B}_{h}^{\mathrm{rot}} = \mathbf{B}_{h}^{n} B h rot = B h n Δ t = 0.025 \Delta t = 0.025 Δ t = 0.025 V = 0 V = 0 V = 0 Δ t = 0.0125 \Delta t = 0.0125 Δ t = 0.0125 B h r o t = B h n \mathbf{B}_{h}^{\mathrm{rot}} = \mathbf{B}_{h}^{n} B h rot = B h n
with m i / m e = 25 m_i/m_e = 25 m i / m e = 25 c 2 = 400 c^2 = 400 c 2 = 400 § 5.4 . The periodic domain of length L = 8 L = 8 L = 8 t = 40 t = 40 t = 40 T = 0.1 T = 0.1 T = 0.1 V = 0.5 V = 0.5 V = 0.5 Δ t = 0.025 \Delta t = 0.025 Δ t = 0.025 0.0125 0.0125 0.0125 Δ t = 0.025 \Delta t = 0.025 Δ t = 0.025 V = 0 V = 0 V = 0
At Δ t = 0.025 \Delta t = 0.025 Δ t = 0.025 B h r o t = B h n \mathbf{B}_{h}^{\mathrm{rot}} = \mathbf{B}_{h}^{n} B h rot = B h n Δ t = 0.0125 \Delta t = 0.0125 Δ t = 0.0125
Consider a uniform, neutral and cold plasma whose species, with the mass densities ρ s \rho_s ρ s κ s \kappa_s κ s V x ^ V\hat{x} V x ^ B 0 z ^ B_0\hat{z} B 0 z ^ V B 0 y ^ VB_0\hat{y} V B 0 y ^ e i k y e^{iky} e ik y ξ s \xi_s ξ s y ^ \hat{y} y ^ u s x ^ + w s y ^ u_s\hat{x} + w_s\hat{y} u s x ^ + w s y ^ E x x ^ + E y y ^ E_x\hat{x} + E_y\hat{y} E x x ^ + E y y ^ b z ^ b\hat{z} b z ^ Φ \boldsymbol{\Phi} Φ S S S
ξ s n + 1 / 2 = ξ s n + δ w s n , ξ s n + 1 = ξ s n + 1 / 2 + δ w s n + 1 , \xi_s^{n+1/2} = \xi_s^n + \delta w_s^n,\qquad\xi_s^{n+1} = \xi_s^{n+1/2} + \delta w_s^{n+1}, ξ s n + 1/2 = ξ s n + δ w s n , ξ s n + 1 = ξ s n + 1/2 + δ w s n + 1 , u ˉ s − u s n = δ κ s ( E ˉ x + B 0 w ˉ s ) , w ˉ s − w s n = δ κ s ( E ˉ y − B 0 u ˉ s − V b r o t ) , \bar{u}_s-u_s^n = \delta\kappa_s\left(\bar{E}_x+B_0\bar{w}_s\right),\qquad\bar{w}_s-w_s^n = \delta\kappa_s\left(\bar{E}_y-B_0\bar{u}_s-Vb^{\mathrm{rot}}\right), u ˉ s − u s n = δ κ s ( E ˉ x + B 0 w ˉ s ) , w ˉ s − w s n = δ κ s ( E ˉ y − B 0 u ˉ s − V b rot ) , ϵ 0 ( E ˉ x − E x n ) = δ [ i k b ˉ μ 0 − ∑ s κ s ρ s ( u ˉ s − i k V ξ s n + 1 / 2 ) ] , \epsilon_0\left(\bar{E}_x-E_x^n\right)=\delta\left[\frac{ik\bar{b}}{\mu_0}-\sum_s\kappa_s\rho_s\left(\bar{u}_s-ikV\xi_s^{n+1/2}\right)\right], ϵ 0 ( E ˉ x − E x n ) = δ [ μ 0 ik b ˉ − s ∑ κ s ρ s ( u ˉ s − ik V ξ s n + 1/2 ) ] , ϵ 0 ( E ˉ y − E y n ) = − δ J ˉ y , b ˉ = b n + δ i k E ˉ x , \epsilon_0\left(\bar{E}_y-E_y^n\right)=-\delta\bar{J}_y,\qquad\bar{b}=b^n+\delta ik\bar{E}_x, ϵ 0 ( E ˉ y − E y n ) = − δ J ˉ y , b ˉ = b n + δ ik E ˉ x , where b r o t = b n + δ i k E x n b^{\mathrm{rot}}=b^n+\delta ikE_x^n b rot = b n + δ ik E x n J y = ∑ s κ s ρ s w s J_y=\sum_s\kappa_s\rho_sw_s J y = ∑ s κ s ρ s w s a ˉ = 1 2 ( a n + a n + 1 ) \bar{a}=\frac{1}{2}(a^n+a^{n+1}) a ˉ = 2 1 ( a n + a n + 1 ) a a a ϵ 0 E y = − ∑ s κ s ρ s ξ s \epsilon_0E_y=-\sum_s\kappa_s\rho_s\xi_s ϵ 0 E y = − ∑ s κ s ρ s ξ s ∑ s κ s ρ s ξ s n + 1 / 2 = − ϵ 0 E y n + δ J y n \sum_s\kappa_s\rho_s\xi_s^{n+1/2}=-\epsilon_0E_y^n+\delta J_y^n ∑ s κ s ρ s ξ s n + 1/2 = − ϵ 0 E y n + δ J y n b r o t − b ˉ = b^{\mathrm{rot}}-\bar{b}= b rot − b ˉ = − 1 2 δ i k ( E x n + 1 − E x n ) -\frac{1}{2}\delta i k(E_x^{n+1}-E_x^n) − 2 1 δ ik ( E x n + 1 − E x n ) E y n − E ˉ y = δ J ˉ y / ϵ 0 E_y^n-\bar{E}_y=\delta\bar{J}_y/\epsilon_0 E y n − E ˉ y = δ J ˉ y / ϵ 0
T ( y n + 1 − y n ) = Δ t 2 L ( y n + 1 + y n ) , y = ( u s , w s , E x , E y , b ) , \mathbf{T}(\mathbf{y}^{n+1}-\mathbf{y}^n)=\frac{\Delta t}{2}\mathbf{L}(\mathbf{y}^{n+1}+\mathbf{y}^n), \qquad\mathbf{y}=(u_s,w_s,E_x,E_y,b), T ( y n + 1 − y n ) = 2 Δ t L ( y n + 1 + y n ) , y = ( u s , w s , E x , E y , b ) , where L L L
d u s d t = κ s ( E x + B 0 w s ) , d w s d t = κ s ( E y − B 0 u s − V b ) , d b d t = i k E x , \frac{\mathrm{d}u_s}{\mathrm{d}t}=\kappa_s(E_x+B_0w_s), \qquad\frac{\mathrm{d}w_s}{\mathrm{d}t}=\kappa_s(E_y-B_0u_s-Vb), \qquad\frac{\mathrm{d}b}{\mathrm{d}t}=ikE_x, d t d u s = κ s ( E x + B 0 w s ) , d t d w s = κ s ( E y − B 0 u s − V b ) , d t d b = ik E x , ϵ 0 d E x d t = i k b μ 0 − ∑ s κ s ρ s u s − i k V ϵ 0 E y , ϵ 0 d E y d t = − J y , \epsilon_0\frac{\mathrm{d}E_x}{\mathrm{d}t}=\frac{ikb}{\mu_0}-\sum_s\kappa_s\rho_su_s-ikV\epsilon_0E_y, \qquad\epsilon_0\frac{\mathrm{d}E_y}{\mathrm{d}t}=-J_y, ϵ 0 d t d E x = μ 0 ik b − s ∑ κ s ρ s u s − ik V ϵ 0 E y , ϵ 0 d t d E y = − J y , and T \mathbf{T} T
( T y ) w s = w s − i k δ 2 V κ s E x , ( T y ) E x = E x + i k δ 2 V ϵ 0 J y . (\mathbf{T}\mathbf{y})_{w_s}=w_s-ik\delta^2V\kappa_sE_x, \qquad(\mathbf{T}\mathbf{y})_{E_x}=E_x+\frac{ik\delta^2V}{\epsilon_0}J_y. ( Ty ) w s = w s − ik δ 2 V κ s E x , ( Ty ) E x = E x + ϵ 0 ik δ 2 V J y . The cold-fluid equations conserve
H 0 ( y ) = ∑ s ρ s ( ∣ u s ∣ 2 + ∣ w s ∣ 2 ) + ϵ 0 ∣ E x ∣ 2 + ϵ 0 ∣ E y − V b ∣ 2 + ( 1 μ 0 − ϵ 0 V 2 ) ∣ b ∣ 2 = y H H 0 y , \mathcal{H}_0(\mathbf{y})=\sum_s\rho_s\left(|u_s|^2+|w_s|^2\right)+\epsilon_0|E_x|^2+\epsilon_0|E_y-Vb|^2+\left(\frac{1}{\mu_0}-\epsilon_0V^2\right)|b|^2=\mathbf{y}^{H}\mathbf{H}_0\mathbf{y}, H 0 ( y ) = s ∑ ρ s ( ∣ u s ∣ 2 + ∣ w s ∣ 2 ) + ϵ 0 ∣ E x ∣ 2 + ϵ 0 ∣ E y − V b ∣ 2 + ( μ 0 1 − ϵ 0 V 2 ) ∣ b ∣ 2 = y H H 0 y , that is, Re ( y H H 0 L y ) = 0 \operatorname{Re}(\mathbf{y}^{H}\mathbf{H}_0\mathbf{L}\mathbf{y})=0 Re ( y H H 0 Ly ) = 0 y \mathbf{y} y T \mathbf{T} T H 0 \mathbf{H}_0 H 0 − i k δ 2 V κ s ρ s -ik\delta^2V\kappa_s\rho_s − ik δ 2 V κ s ρ s w s w_s w s i k δ 2 V κ s ρ s ik\delta^2V\kappa_s\rho_s ik δ 2 V κ s ρ s E x E_x E x H 0 T \mathbf{H}_0\mathbf{T} H 0 T H Δ t ( y ) = y H H 0 T y \mathcal{H}_{\Delta t}(\mathbf{y})=\mathbf{y}^{H}\mathbf{H}_0\mathbf{T}\mathbf{y} H Δ t ( y ) = y H H 0 Ty
H Δ t ( y n + 1 ) − H Δ t ( y n ) = Re [ ( y n + 1 + y n ) H H 0 T ( y n + 1 − y n ) ] = Δ t 2 Re [ ( y n + 1 + y n ) H H 0 L ( y n + 1 + y n ) ] = 0. \begin{aligned}
\mathcal{H}_{\Delta t}(\mathbf{y}^{n+1})-\mathcal{H}_{\Delta t}(\mathbf{y}^{n})=\operatorname{Re}\left[(\mathbf{y}^{n+1}+\mathbf{y}^{n})^{H}\mathbf{H}_0\mathbf{T}(\mathbf{y}^{n+1}-\mathbf{y}^{n})\right] \\
&=\frac{\Delta t}{2}\operatorname{Re}\left[(\mathbf{y}^{n+1}+\mathbf{y}^{n})^{H}\mathbf{H}_0\mathbf{L}(\mathbf{y}^{n+1}+\mathbf{y}^{n})\right]=0.
\end{aligned} H Δ t ( y n + 1 ) − H Δ t ( y n ) = Re [ ( y n + 1 + y n ) H H 0 T ( y n + 1 − y n ) ] = 2 Δ t Re [ ( y n + 1 + y n ) H H 0 L ( y n + 1 + y n ) ] = 0. Completing the square gives
H Δ t ( y ) = ∑ s ρ s ∣ u s ∣ 2 + ∑ s ρ s ∣ w s − i k δ 2 V κ s E x ∣ 2 + ϵ 0 [ 1 − ( k V ω p δ 2 ) 2 ] ∣ E x ∣ 2 + ϵ 0 ∣ E y − V b ∣ 2 + ( 1 μ 0 − ϵ 0 V 2 ) ∣ b ∣ 2 , \begin{aligned}
\mathcal{H}_{\Delta t}(\mathbf{y})=\sum_s\rho_s|u_s|^2+\sum_s\rho_s\left|w_s-ik\delta^2V\kappa_sE_x\right|^2+\epsilon_0\left[1-(kV\omega_p\delta^2)^2\right]|E_x|^2 \\
&+\epsilon_0|E_y-Vb|^2+\left(\frac{1}{\mu_0}-\epsilon_0V^2\right)|b|^2,
\end{aligned} H Δ t ( y ) = s ∑ ρ s ∣ u s ∣ 2 + s ∑ ρ s w s − ik δ 2 V κ s E x 2 + ϵ 0 [ 1 − ( k V ω p δ 2 ) 2 ] ∣ E x ∣ 2 + ϵ 0 ∣ E y − V b ∣ 2 + ( μ 0 1 − ϵ 0 V 2 ) ∣ b ∣ 2 , with ω p 2 = ϵ 0 − 1 ∑ s κ s 2 ρ s \omega_p^2=\epsilon_0^{-1}\sum_s\kappa_s^2\rho_s ω p 2 = ϵ 0 − 1 ∑ s κ s 2 ρ s ∣ V ∣ < c |V|<c ∣ V ∣ < c H Δ t \mathcal{H}_{\Delta t} H Δ t ∣ k V ∣ ω p Δ t 2 < 4 |kV|\omega_p\Delta t^2<4 ∣ k V ∣ ω p Δ t 2 < 4 ∣ k ∣ ≤ k ∗ |k|\le k_* ∣ k ∣ ≤ k ∗ k ∗ ∣ V ∣ ω p Δ t 2 < 4 k_*|V|\omega_p\Delta t^2<4 k ∗ ∣ V ∣ ω p Δ t 2 < 4 β i < 4 \beta_i<4 β i < 4 ω p e Δ t ≈ 2.5 \omega_{pe}\Delta t\approx2.5 ω p e Δ t ≈ 2.5 1.25 1.25 1.25 β i ≈ 0.76 \beta_i\approx0.76 β i ≈ 0.76 0.19 0.19 0.19
The entry of T \mathbf{T} T E x E_x E x w s w_s w s b r o t = b n b^{\mathrm{rot}}=b^n b rot = b n w s w_s w s w s n + 1 − w s n = Δ t ( L y ˉ ) w s + 2 i k δ 2 V κ s E ˉ x w_s^{n+1}-w_s^n=\Delta t(L\bar{\mathbf{y}})_{w_s}+2ik\delta^2V\kappa_s\bar{E}_x w s n + 1 − w s n = Δ t ( L y ˉ ) w s + 2 ik δ 2 V κ s E ˉ x H 0 T \mathbf{H}_0\mathbf{T} H 0 T H 0 \mathcal{H}_0 H 0 H Δ t \mathcal{H}_{\Delta t} H Δ t Δ t = 0.025 \Delta t=0.025 Δ t = 0.025
Proof of Proposition 2
Proof. (i) For a charged species s s s Ω i \Omega_i Ω i ζ s , i \zeta_{s,i} ζ s , i ζ s , i , ℓ \zeta_{s,i,\ell} ζ s , i , ℓ Ω i \Omega_i Ω i ∑ s , i k B T s , i 2 m s ζ s , i T Γ s , i ζ s , i \sum_{s,i}\frac{k_B T_{s,i}}{2m_s}\zeta_{s,i}^{\mathsf T}\Gamma_{s,i}\zeta_{s,i} ∑ s , i 2 m s k B T s , i ζ s , i T Γ s , i ζ s , i
( Γ s , i ) ℓ ℓ ′ = ∑ k ∈ P s , i m k λ ℓ ( x k ) λ ℓ ′ ( x k ) + ∫ Ω i ρ s h , c λ ℓ λ ℓ ′ d x , (\Gamma_{s,i})_{\ell\ell'}=\sum_{k\in\mathcal{P}_{s,i}}m_k\lambda_\ell(\boldsymbol{x}_k)\lambda_{\ell'}(\boldsymbol{x}_k)+\int_{\Omega_i}\rho_s^{h,c}\lambda_\ell\lambda_{\ell'}\,\mathrm{d}x, ( Γ s , i ) ℓ ℓ ′ = k ∈ P s , i ∑ m k λ ℓ ( x k ) λ ℓ ′ ( x k ) + ∫ Ω i ρ s h , c λ ℓ λ ℓ ′ d x , the Gram matrix of the nodal functions λ ℓ \lambda_\ell λ ℓ Ω i \Omega_i Ω i
The problem has a feasible point. As r ∥ \boldsymbol{r}_{\parallel} r ∥ D D D e F \boldsymbol{e}_F e F ϵ 0 D e F = ϵ 0 D e c − r ∥ \epsilon_0D\boldsymbol{e}_F=\epsilon_0D\boldsymbol{e}^c-\boldsymbol{r}_{\parallel} ϵ 0 D e F = ϵ 0 D e c − r ∥ ζ = 0 \zeta=0 ζ = 0 Q n + 1 = Q \boldsymbol{Q}^{n+1}=\boldsymbol{Q} Q n + 1 = Q W s , i n + 1 = W s , i c \boldsymbol{W}_{s,i}^{n+1}=\boldsymbol{W}_{s,i}^c W s , i n + 1 = W s , i c ϵ 0 D e F − Q n + 1 = r − r ∥ = r ⊥ \epsilon_0D\boldsymbol{e}_F-\boldsymbol{Q}^{n+1}=\boldsymbol{r}-\boldsymbol{r}_{\parallel}=\boldsymbol{r}_{\perp} ϵ 0 D e F − Q n + 1 = r − r ∥ = r ⊥ ( e F , 0 ) (\boldsymbol{e}_F,0) ( e F , 0 )
The cost depends on ζ \zeta ζ η s , i = Γ s , i 1 / 2 ζ s , i \eta_{s,i}=\Gamma_{s,i}^{1/2}\zeta_{s,i} η s , i = Γ s , i 1/2 ζ s , i ( e n + 1 , η ) (\boldsymbol{e}^{n+1},\eta) ( e n + 1 , η )
J = ϵ 0 2 ( e n + 1 − e c ) T M E ( e n + 1 − e c ) + ∑ s , i k B T s , i 2 m s ∣ η s , i ∣ 2 , \mathcal{J}=\frac{\epsilon_0}{2}(\boldsymbol{e}^{n+1}-\boldsymbol{e}^c)^{\mathsf T}M_E(\boldsymbol{e}^{n+1}-\boldsymbol{e}^c)+\sum_{s,i}\frac{k_B T_{s,i}}{2m_s}|\eta_{s,i}|^2, J = 2 ϵ 0 ( e n + 1 − e c ) T M E ( e n + 1 − e c ) + s , i ∑ 2 m s k B T s , i ∣ η s , i ∣ 2 , which is strictly convex and coercive, as M E M_E M E T s , i > 0 T_{s,i}>0 T s , i > 0 ( e n + 1 , ζ ) ↦ ( e n + 1 , η ) (\boldsymbol{e}^{n+1},\zeta)\mapsto(\boldsymbol{e}^{n+1},\eta) ( e n + 1 , ζ ) ↦ ( e n + 1 , η ) J \mathcal{J} J e n + 1 \boldsymbol{e}^{n+1} e n + 1 η \eta η ζ s , i \zeta_{s,i} ζ s , i y y y Γ s , i y = 0 \Gamma_{s,i}y=0 Γ s , i y = 0
This difference does not change the corrected state. The state (27) and the charges Q n + 1 \boldsymbol{Q}^{n+1} Q n + 1 ζ s , i \zeta_{s,i} ζ s , i ζ s \zeta_s ζ s P s , i \mathcal{P}_{s,i} P s , i ζ s h s c \zeta_s h_s^c ζ s h s c Ω i \Omega_i Ω i z = ∑ ℓ y ℓ λ ℓ z=\sum_\ell y_\ell\lambda_\ell z = ∑ ℓ y ℓ λ ℓ
0 = y T Γ s , i y = ∑ k ∈ P s , i m k z ( x k ) 2 + ∫ Ω i ρ s h , c z 2 d x , 0=y^{\mathsf T}\Gamma_{s,i}y=\sum_{k\in\mathcal{P}_{s,i}}m_k z(\boldsymbol{x}_k)^2+\int_{\Omega_i}\rho_s^{h,c}z^2\,\mathrm{d}x, 0 = y T Γ s , i y = k ∈ P s , i ∑ m k z ( x k ) 2 + ∫ Ω i ρ s h , c z 2 d x , so that z z z P s , i \mathcal{P}_{s,i} P s , i ρ s h , c > 0 \rho_s^{h,c}>0 ρ s h , c > 0 h s c = M [ W s h , c ] h_s^c=\mathcal{M}[\boldsymbol{W}_s^{h,c}] h s c = M [ W s h , c ] ρ s h , c \rho_s^{h,c} ρ s h , c z h s c zh_s^c z h s c Ω i \Omega_i Ω i
(ii) Every solution satisfies (29). Its first equation states that the Gauss defect after the correction is r ⊥ \boldsymbol{r}_{\perp} r ⊥
(iii) The total energy is the sum of the energies of the species and the electric and magnetic energies. The correction leaves b \boldsymbol{b} b
