Structure-preserving strategy for conservative simulation of relativistic nonlinear Landau--Fokker--Planck equation

TL;DR

This paper introduces a structure-preserving numerical strategy for the relativistic Landau--Fokker--Planck equation that maintains its fundamental symmetries, enabling accurate mass, momentum, and energy conservation in simulations.

Contribution

The paper presents a novel simulation method that preserves the intrinsic symmetries of the Beliaev--Budker kernel, ensuring physical conservation laws without artificial constraints.

Findings

Conservation of mass, momentum, and energy achieved in simulations.

Symmetry preservation improves the physical accuracy of relativistic plasma modeling.

Method demonstrates robustness and fidelity in relativistic nonlinear Landau--Fokker--Planck simulations.

Abstract

Mathematical symmetries of the Beliaev--Budker kernel are the most important structure of the relativistic Landau--Fokker--Planck equation. By preserving the beautiful symmetries, a mass-momentum-energy-conserving simulation has been demonstrated without any artificial constraints.

Figures

Equations

$$\frac{\partial f_{\mathrm{s}}}{\partial t}=\frac{\Gamma }{2m_{\mathrm{s}}^2}\frac{\partial }{\partial \mathbf{u}}\cdot \int \mathsf{U}(\mathbf{u},{\mathbf{u}}^{\prime })\cdot \left(f_{{\mathrm{s}}^{\prime }}\frac{\partial f_{\mathrm{s}}}{\partial \mathbf{u}}-\frac{m_{\mathrm{s}}}{m_{{\mathrm{s}}^{\prime }}}{f}_{\mathrm{s}}\frac{\partial f_{{\mathrm{s}}^{\prime }}}{\partial {\mathbf{u}}^{\prime }}\right)\mathrm{d}{\mathbf{u}}^{\prime },$$

$$\mathsf{U}=\frac{r^2}{\gamma {\gamma }^{\prime }{w}^3}\left[w^2\mathsf{I}-{\gamma }^2\mathbf{vv}-{\gamma }^{\prime 2}{\mathbf{v}}^{\prime }{\mathbf{v}}^{\prime }+r\gamma {\gamma }^{\prime }(\mathbf{v}{\mathbf{v}}^{\prime }+{\mathbf{v}}^{\prime }\mathbf{v})\right],$$

$$\gamma =\sqrt{1+|\mathbf{u}/c{|}^2}=1/\sqrt{1-|\mathbf{v}/c{|}^2},$$

$${\gamma }^{\prime }=\sqrt{1+|{\mathbf{u}}^{\prime }/c{|}^2}=1/\sqrt{1-|{\mathbf{v}}^{\prime }/c{|}^2}.$$

$$f_{\mathrm{s}}{f}_{{\mathrm{s}}^{\prime }}\mathsf{U}(\mathbf{u},{\mathbf{u}}^{\prime })\cdot \left(\frac{\partial \log f_{\mathrm{s}}}{\partial \mathbf{u}}-\frac{m_{\mathrm{s}}}{m_{{\mathrm{s}}^{\prime }}}\frac{\partial \log f_{{\mathrm{s}}^{\prime }}}{\partial {\mathbf{u}}^{\prime }}\right).$$

$$\frac{\partial {f_{\mathrm{s}}}^{\mathbf{j}}}{\partial t}\equiv \frac{\Gamma }{2{m_{\mathrm{s}}}^2}\frac{\delta }{\delta \mathbf{u}}\cdot \sum _{\mathbf{k}}$$

$$\left({f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}\frac{\delta {f_{\mathrm{s}}}^{\mathbf{j}}}{\delta \mathbf{u}}-\frac{m_{\mathrm{s}}}{m_{{\mathrm{s}}^{\prime }}}{f_{\mathrm{s}}}^{\mathbf{j}}\frac{\delta {f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}\right)\Delta {\mathbf{u}}^{\prime },$$

$$\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta u_1}\equiv \frac{a^{j_1+1,j_2,j_3,k_1,k_2,k_3}-a^{j_1-1,j_2,j_3,k_1,k_2,k_3}}{2\Delta u_1},$$

$$\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta u_2}\equiv \frac{a^{j_1,j_2+1,j_3,k_1,k_2,k_3}-a^{j_1,j_2-1,j_3,k_1,k_2,k_3}}{2\Delta u_2},$$

$$\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta u_3}\equiv \frac{a^{j_1,j_2,j_3+1,k_1,k_2,k_3}-a^{j_1,j_2,j_3-1,k_1,k_2,k_3}}{2\Delta u_3},$$

$$\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta u_1^{\prime }}\equiv \frac{a^{j_1,j_2,j_3,k_1+1,k_2,k_3}-a^{j_1,j_2,j_3,k_1-1,k_2,k_3}}{2\Delta u_1^{\prime }},$$

$$\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta u_2^{\prime }}\equiv \frac{a^{j_1,j_2,j_3,k_1,k_2+1,k_3}-a^{j_1,j_2,j_3,k_1,k_2-1,k_3}}{2\Delta u_2^{\prime }},$$

$$\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta u_3^{\prime }}\equiv \frac{a^{j_1,j_2,j_3,k_1,k_2,k_3+1}-a^{j_1,j_2,j_3,k_1,k_2,k_3-1}}{2\Delta u_3^{\prime }}.$$

$$\mathsf{U}({\mathbf{u}}^{\mathbf{j}},{{\mathbf{u}}^{\prime }}^{\mathbf{k}})\equiv \frac{(r^{\mathbf{j},\mathbf{k}}{)}^2}{{\gamma }^{\mathbf{j}}{{\gamma }^{\prime }}^{\mathbf{k}}(w^{\mathbf{j},\mathbf{k}}{)}^3}[(w^{\mathbf{j},\mathbf{k}}{)}^2\mathsf{I}-{\mathbf{u}}^{\mathbf{j}}{\mathbf{u}}^{\mathbf{j}}-{{\mathbf{u}}^{\prime }}^{\mathbf{k}}{{\mathbf{u}}^{\prime }}^{\mathbf{k}}$$

$$+r^{\mathbf{j},\mathbf{k}}({\mathbf{u}}^{\mathbf{j}}{{\mathbf{u}}^{\prime }}^{\mathbf{k}}+{{\mathbf{u}}^{\prime }}^{\mathbf{k}}{\mathbf{u}}^{\mathbf{j}})],$$

$$\mathsf{U}({\mathbf{u}}^{\mathbf{j}},{{\mathbf{u}}^{\prime }}^{\mathbf{k}})=\mathsf{U}({{\mathbf{u}}^{\prime }}^{\mathbf{k}},{\mathbf{u}}^{\mathbf{j}}),$$

$$\mathsf{U}({\mathbf{u}}^{\mathbf{j}},{{\mathbf{u}}^{\prime }}^{\mathbf{k}})\cdot {\mathbf{v}}^{\mathbf{j}}=\mathsf{U}({{\mathbf{u}}^{\prime }}^{\mathbf{k}},{\mathbf{u}}^{\mathbf{j}})\cdot {{\mathbf{v}}^{\prime }}^{\mathbf{k}}.$$

$$\sum _{\mathbf{j},\mathbf{k}}\left(\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta \mathbf{u}}{b}^{\mathbf{j},\mathbf{k}}\right)\Delta \mathbf{u}\Delta {\mathbf{u}}^{\prime }=-\sum _{\mathbf{j},\mathbf{k}}\left(a^{\mathbf{j},\mathbf{k}}\frac{\delta b^{\mathbf{j},\mathbf{k}}}{\delta \mathbf{u}}\right)\Delta \mathbf{u}\Delta {\mathbf{u}}^{\prime },$$

$$\sum _{\mathbf{j},\mathbf{k}}\left(\frac{\delta a^{\mathbf{j},\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}{b}^{\mathbf{j},\mathbf{k}}\right)\Delta \mathbf{u}\Delta {\mathbf{u}}^{\prime }=-\sum _{\mathbf{j},\mathbf{k}}\left(a^{\mathbf{j},\mathbf{k}}\frac{\delta b^{\mathbf{j},\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}\right)\Delta \mathbf{u}\Delta {\mathbf{u}}^{\prime },$$

$$\frac{\partial }{\partial t}\left(\sum _{\mathbf{j}}{m}_{\mathrm{s}}{\mathbf{u}}^{\mathbf{j}}{f_{\mathrm{s}}}^{\mathbf{j}}\Delta \mathbf{u}\right)$$

$$=$$

$$\left({f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}\frac{\delta {f_{\mathrm{s}}}^{\mathbf{j}}}{\delta \mathbf{u}}-\frac{m_{\mathrm{s}}}{m_{{\mathrm{s}}^{\prime }}}{f_{\mathrm{s}}}^{\mathbf{j}}\frac{\delta {f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}\right)\Delta {\mathbf{u}}^{\prime }]\Delta \mathbf{u}$$

$$=$$

$$\left(\because \delta {\mathbf{u}}^{\mathbf{j}}/\delta \mathbf{u}\text{mbox}isanidentitymatrix.\right)$$

$$\frac{\partial }{\partial t}\left(\sum _{\mathbf{k}}{m}_{{\mathrm{s}}^{\prime }}{{\mathbf{u}}^{\prime }}^{\mathbf{k}}{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}\Delta {\mathbf{u}}^{\prime }\right)$$

$$=$$

$$\frac{\partial }{\partial t}\left(\sum _{\mathbf{j}}{m}_{\mathrm{s}}{\mathbf{u}}^{\mathbf{j}}{f_{\mathrm{s}}}^{\mathbf{j}}\Delta \mathbf{u}+\sum _{\mathbf{k}}{m}_{{\mathrm{s}}^{\prime }}{{\mathbf{u}}^{\prime }}^{\mathbf{k}}{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}\Delta {\mathbf{u}}^{\prime }\right)$$

$$=$$

$$\left(\frac{{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{m_{\mathrm{s}}}\frac{\delta {f_{\mathrm{s}}}^{\mathbf{j}}}{\delta \mathbf{u}}-\frac{{f_{\mathrm{s}}}^{\mathbf{j}}}{m_{{\mathrm{s}}^{\prime }}}\frac{\delta {f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}\right)]\Delta \mathbf{u}\Delta {\mathbf{u}}^{\prime }.$$

$$\frac{\partial }{\partial t}\left(\sum _{\mathbf{j}}{m}_{\mathrm{s}}{c}^2{\gamma }^{\mathbf{j}}{f_{\mathrm{s}}}^{\mathbf{j}}\Delta \mathbf{u}\right)$$

$$=$$

$$\left({f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}\frac{\delta {f_{\mathrm{s}}}^{\mathbf{j}}}{\delta \mathbf{u}}-\frac{m_{\mathrm{s}}}{m_{{\mathrm{s}}^{\prime }}}{f_{\mathrm{s}}}^{\mathbf{j}}\frac{\delta {f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}\right)\Delta {\mathbf{u}}^{\prime }]\Delta \mathbf{u}$$

$$=$$

$$\left(\frac{{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{m_{\mathrm{s}}}\frac{\delta {f_{\mathrm{s}}}^{\mathbf{j}}}{\delta \mathbf{u}}-\frac{{f_{\mathrm{s}}}^{\mathbf{j}}}{m_{{\mathrm{s}}^{\prime }}}\frac{\delta {f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}\right)]\Delta \mathbf{u}\Delta {\mathbf{u}}^{\prime }.$$

$$\frac{\partial }{\partial t}\left(\sum _{\mathbf{k}}{m}_{{\mathrm{s}}^{\prime }}{c}^2{{\gamma }^{\prime }}^{\mathbf{k}}{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}\Delta {\mathbf{u}}^{\prime }\right)$$

$$=$$

$$\left(\frac{{f_{\mathrm{s}}}^{\mathbf{j}}}{m_{{\mathrm{s}}^{\prime }}}\frac{\delta {f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}-\frac{{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{m_{\mathrm{s}}}\frac{\delta {f_{\mathrm{s}}}^{\mathbf{j}}}{\delta \mathbf{u}}\right)]\Delta {\mathbf{u}}^{\prime }\Delta \mathbf{u}.$$

$$\frac{\partial }{\partial t}\left(\sum _{\mathbf{j}}{m}_{\mathrm{s}}{c}^2{\gamma }^{\mathbf{j}}{f_{\mathrm{s}}}^{\mathbf{j}}\Delta \mathbf{u}+\sum _{\mathbf{k}}{m}_{{\mathrm{s}}^{\prime }}{c}^2{{\gamma }^{\prime }}^{\mathbf{k}}{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}\Delta {\mathbf{u}}^{\prime }\right)$$

$$=$$

$$\left(\frac{{f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{m_{\mathrm{s}}}\frac{\delta {f_{\mathrm{s}}}^{\mathbf{j}}}{\delta \mathbf{u}}-\frac{{f_{\mathrm{s}}}^{\mathbf{j}}}{m_{{\mathrm{s}}^{\prime }}}\frac{\delta {f_{{\mathrm{s}}^{\prime }}}^{\mathbf{k}}}{\delta {\mathbf{u}}^{\prime }}\right)]\Delta \mathbf{u}\Delta {\mathbf{u}}^{\prime },$$

$$\bar{\mathbf{v}}\equiv \frac{\delta ({\gamma }^{\mathbf{j}}{c}^2)}{\delta \mathbf{u}},{\bar{\mathbf{v}}}^{\prime }\equiv \frac{\delta ({{\gamma }^{\prime }}^{\mathbf{k}}{c}^2)}{\delta {\mathbf{u}}^{\prime }}.$$

$$({\gamma }^{j_1+1,j_2,j_3}{)}^2-({\gamma }^{j_1-1,j_2,j_3}{)}^2=({u_1}^{j_1+1}/c{)}^2-({u_1}^{j_1-1}/c{)}^2,$$

Full text

I Introduction

The Landau–Fokker–Planck (LFP) equation in the International System of Units is Lifshitz and Pitaevskii (1981)

[TABLE]

where $\Gamma=(q_{\mathrm{s}}^{2}q_{\mathrm{s^{\prime}}}^{2}\log\Lambda)/4\pi{\varepsilon_{0}}^{2}$, $\log\Lambda$ is the Coulomb logarithm, $\varepsilon_{0}$ is the vacuum permittivity, $(f_{\mathrm{s}},f_{\mathrm{s^{\prime}}})$, $(m_{\mathrm{s}},m_{\mathrm{s^{\prime}}})$, $(q_{\mathrm{s}},q_{\mathrm{s^{\prime}}})$ and $(\mathbf{u},\mathbf{u^{\prime}})$ are the distribution function, mass, electric charge and momentum per unit mass of species $(\mathrm{s},\mathrm{s^{\prime}})$, respectively. In the relativistic case, the collision kernel $\mathsf{U}(\mathbf{u},\mathbf{u^{\prime}})$ is described as follows Beliaev and Budker (1956):

[TABLE]

where $r=\gamma\gamma^{\prime}(1-\mathbf{v}\cdot\mathbf{v^{\prime}}/c^{2})$, $w=c\sqrt{r^{2}-1}$, $\mathbf{v}=\mathbf{u}/\gamma$, $\mathbf{v^{\prime}}=\mathbf{u^{\prime}}/\gamma^{\prime}$, and $c$ is the speed of light in vacuum. The Lorentz factor is defined as follows:

[TABLE]

The relativistic LFP equation is designed so as to ensure the mass-momentum-energy conservation and the H-theorem.

The collisional processes described by the relativistic LFP equation are essential in fusion plasmas. In practical simulations, the relativistic LFP equation is often linearized so that the computational cost is reduced from $O(N^{2})$ to $O(N)$, where $N$ is the number of unknowns. For fast ignition, which heats the imploded core by relativistic fast electrons, of inertial confinement fusion Tabak et al. (1994), a linearization of Nakashima and Takabe Nakashima and Takabe (2002) is employed by relativistic Fokker–Planck codes such as RFP-2D Yokota et al. (2006), FIBMET Johzaki et al. (2009) and FIDO Sherlock (2009); Ridgers et al. (2011). The linearization is based on the fact that the colliding particles are much faster than the collided ones. This violates symmetry of the collision kernel so the conservation laws are maintained only at the continuous limit. For runaway electrons in tokamak disruptions Dreicer (1959), linearization which assumes an weakly relativistic equilibrium background Karney and Fisch (1985); Karney (1986) is sometimes performed to take into account the effect of non-thermal electrons Nuga et al. (2016), while the conservation laws are violated. Further, TASK/FP Nuga and Fukuyama (2011) and CQL3D Petrov and Harvey (2016) codes have options that decompose the nonlinear LFP equation into Legendre modes and solve first a few modes. The Legendre polynomials ensure the mass-momentum-energy conservation limit to resolve structures of the pitch angle. Recently, Stahl et al. developed NORSE code Stahl et al. (2017) which models nonlinear electron–electron collisions by the Braams–Karney potential formulation Braams and Karnel (1987). However, the NORSE code violates the conservation laws including the mass conservation. In the potential form, the “nonlinear constraints” is one of the ways to preserve the conservation laws Taitano et al. (2015) and equilibrium state Taitano et al. (2017), i.e. one of the projection method in the literature of applied mathematics Eich (1993). However, the projection method can affect the stability of numerical schemes, and it may not be suit for long time-scale simulations. In magnetohydrodynamics simulations, for example, the projection method Brackbill and Barnes (1980) or constrained transport (CT) method Evans and Hawley (1988) has been used to enforce the solenoidal constraint ($\nabla\cdot\mathbf{B}=0$). However, such an inconsistent magnetic field can induce “checkerboard phenomenon” which is one of the numerical instabilities Flock et al. (2010).

Unlike the relativistic regime, many non-relativistic structure-preserving schemes have been proposed to ensure the conservation laws, positivity and H-theorem. Chang and Cooper developed a positivity-preserving scheme for one-dimensional (1D) linearized LFP equation Chang and Cooper (1970), and it was extended to nonlinear isotropic LFP equation Buet and Cordier (2002) and nonlinear multi-dimensional one Yoon and Chang (2014a, b); Hager et al. (2016) preserving the conservation laws and H-theorem. A structure-preserving finite-element scheme Hirvijoki and Adams (2017) is also developed to ensure the conservation laws on unstructure meshes. These works are based on a weak-form associated with Eq. (1). The integrand of Eq. (1) can be transformed into the following one analytically:

[TABLE]

The weak-form coming from Eq. (5) is called as “log” weak-form, and the proof of H-theorem by the log weak-form is straightforward. An analytical discussion was first given by Pekker and Khudik Pekker and Khudik (1984), and conservative and entropic discretizations have been proposed for isotropic Berezin et al. (1987); Buet and Cordier (1998) and multi-dimensional cases Degond and Lucquin-Desreux (1994); Buet et al. (1997). Furthermore, an energy-conserving LFP scheme with the Rosenbluth potential form Rosenbluth et al. (1957) was proposed using an analogy of the Maxwell stress tensor in electromagnetism Chacón et al. (2000a, b). However, structure-preserving schemes for the relativistic LFP equation have not been developed since the Lorentz factor is not expressed as polynomials of finite-order in the momentum Hirvijoki and Adams (2017).

In this paper, we demonstrate a structure-preserving simulation of the relativistic LFP equation which strictly preserve the conservation laws of mass, momentum and energy. The key concept is the same with our recent work about a quadratic conservative scheme for the relativistic Vlasov–Maxwell system Shiroto et al. (2019). Moreover, a similar approach has been used in a non-relativistic scheme with non-uniform meshes Buet and Thanh (2006, 2007). The rest of this paper is composed as follows. Section II shows an intuitive discretization which cannot maintain the energy conservation. Section III deductively derives requirements for the conservation laws in discrete form. Our structure-preserving scheme and its concept are introduced in Sec. IV. A verification of conservation property through a thermal-equilibration is performed in Sec. V. Section VI is conclusions of this article.

II CONVENTIONAL SCHEME

In this article, the relativistic LFP equation is discretized as follows. Note that the conservation laws are not affected by the temporal structure, so only momentum dimensions are discretized here:

[TABLE]

where $\mathbf{j}\equiv[j_{1},j_{2},j_{3}]$ and $\mathbf{k}\equiv[k_{1},k_{2},k_{3}]$ are the indices of uniform momentum grids for $\mathbf{u}$ and $\mathbf{u^{\prime}}$, respectively. $\Delta\mathbf{u}\equiv\Delta u_{1}\Delta u_{2}\Delta u_{3}$, and $(\Delta u_{1},\Delta u_{2},\Delta u_{3})$ is the grid interval of $(u_{1},u_{2},u_{3})$. The species $(\mathrm{s},\mathrm{s^{\prime}})$ use the same momentum grids $(\mathbf{u^{j}}=\mathbf{{u^{\prime}}^{j}},\Delta\mathbf{u}=\Delta\mathbf{u^{\prime}})$. Here we use a second-order central difference for simplicity:

[TABLE]

[TABLE]

where $a^{\mathbf{j},\mathbf{k}}\equiv a(\mathbf{u^{j}},\mathbf{u^{\prime k}})$ is an arbitrary function.

In an intuitive discretization, the collision kernel is calculated by its arguments directly:

[TABLE]

From the definition, the collision kernel Eq. (7) satisfies two mathematical symmetries:

[TABLE]

III REQUIREMENTS FOR CONSERVATION LAWS IN DISCRETE FORM

In this section, requirements for the conservation laws are derived deductively. The mass conservation is trivially maintained, so the discussion is focused on the momentum-energy conservation. The following points are required to prove the conservation laws analytically from the relativistic LFP equation:

Integration-by-parts must be maintained. 2. 2.

$\mathsf{U}(\mathbf{u},\mathbf{u^{\prime}})=\mathsf{U}(\mathbf{u^{\prime}},\mathbf{u})$ is required for the momentum conservation. 3. 3.

$\mathsf{U}(\mathbf{u},\mathbf{u^{\prime}})\cdot\mathbf{v}=\mathsf{U}(\mathbf{u^{\prime}},\mathbf{u})\cdot\mathbf{v^{\prime}}$ is required for the energy conservation.

If one can assume $a^{\mathbf{j,k}}=0$ at the momentum boundaries, the following summation-by-parts, i.e. the integration-by-parts in discrete form, is valid:

[TABLE]

where $b^{\mathbf{j,k}}$ is also an arbitrary function. Therefore, the first point is automatically preserved if a finite-difference operator has a linearity. In addition, the second point is naturally satisfied unless the Braams and Karney potential is employed. The most important discussion in our article is the third point. To satisfy $a^{\mathbf{j,k}}=0$ in the relativistic LFP equation, the computational domain must be large enough to ensure that the distribution function is negligible small at the boundary. In other words, no boundary condition can enforce all of the conservation laws.

III.A MOMENTUM CONSERVATION

The momentum of species “$\mathrm{s}$” is described as a first-order moment of Eq. (6):

[TABLE]

Likewise, the momentum of species “$\mathrm{s^{\prime}}$” is obtained as follows:

[TABLE]

The temporal development of total momentum is obtained as a sum of Eqs. (10) and (11):

[TABLE]

Therefore, the intuitive kernel of Eq. (7) strictly maintains the momentum conservation owing to Eq. (8).

III.B ENERGY CONSERVATION

The energy of species “$\mathrm{s}$” is described as a second-order moment of Eq. (6):

[TABLE]

Likewise, the energy of species “$\mathrm{s^{\prime}}$” is obtained as follows:

[TABLE]

The temporal development of total energy is obtained as a sum of Eqs. (13) and (14):

[TABLE]

where the velocities with overlines are defined as follows:

[TABLE]

Therefore, $\mathsf{U}(\mathbf{u^{j}},\mathbf{{u^{\prime}}^{k}})\cdot\mathbf{\bar{v}}=\mathsf{U}(\mathbf{{u^{\prime}}^{k}},\mathbf{u^{j}})\cdot\mathbf{\bar{v}^{\prime}}$ is required for the energy conservation in discrete form. The intuitive kernel Eq. (7) satisfies Eq. (9), so the energy conservation would be preserved even in discrete form if $\mathbf{\bar{v}}=\mathbf{v^{j}}$ and $\mathbf{\bar{v}^{\prime}}=\mathbf{{v^{\prime}}^{k}}$ were true. Generally speaking, the proposition is false resulting in a violation of the energy conservation. For the second-order central difference;

[TABLE]

The only exception is the non-relativistic limit when the Lorentz factor is always unity. Hirvijoki and Adams reported that their scheme cannot maintain the exact energy conservation in the relativistic regime Hirvijoki and Adams (2017), and our discussion should be connected on the fundamental level with this issue although their discussion was based on the finite-element method.

IV STRUCTURE-PRESERVING STRATEGY

Here we propose a structure-preserving scheme for the relativistic LFP equation which resolves the energy-conservation problem of the intuitive scheme. According to Eqs. (12) and (15), the following discrete requirements must be preserved to maintain the mass-momentum-energy conservation:

Summation-by-parts 2. 2.

$\mathsf{U}(\mathbf{u^{j}},\mathbf{{u^{\prime}}^{k}})=\mathsf{U}(\mathbf{{u^{\prime}}^{k}},\mathbf{u^{j}})$ 3. 3.

$\mathsf{U}(\mathbf{u^{j}},\mathbf{{u^{\prime}}^{k}})\cdot\mathbf{\bar{v}}=\mathsf{U}(\mathbf{{u^{\prime}}^{k}},\mathbf{u^{j}})\cdot\mathbf{\bar{v}^{\prime}}$

The following is the only Beliaev–Budker kernel that preserves the above requirements:

[TABLE]

where the variables with overlines are

[TABLE]

A combination of Eqs. (6) and (17) naturally preserves the law of energy conservation naturally. However, the positivity of distribution function and H-theorem are not guaranteed unconditionally by this formulation. We do not give a discussion of the temporal discretization in this article, which is done in the papers about entropic schemes.

V DEMONSTRATION

To verify the proposed scheme, a collisional relaxation of a particle–antiparticle plasma is calculated for simplicity ($m_{\mathrm{s}}=m_{\mathrm{s^{\prime}}}\equiv m,q_{\mathrm{s}}=-q_{\mathrm{s^{\prime}}}\equiv q$). The particles and antiparticles are initialized by the following shifted Maxwell–Jüttner distrubtions Zenitani (2015).

[TABLE]

where $\theta=0.01$ is the temperature normalized by the rest mass energy, $\mathbf{v}_{0}/c=^{\mathrm{T}}[0.2,0.2,0.2]$, and $\gamma_{0}=1/\sqrt{1-\mathbf{v}_{0}\cdot\mathbf{v}_{0}/c^{2}}$. The computational domain is set to be $\{(u_{1},u_{2},u_{3})|-1.2c\leq u_{1},u_{2},u_{3}\leq 1.8c\}$, and the number of computational cells is $128\times 128\times 128$. In this verification, a first-order Euler explicit method is used as a time integration. The temporal interval is given as $\Gamma m^{-2}\Delta t=1/80$, so that the temporal resolution is fine enough to see the relaxation. In this test, only unlike-particle collisions are considered.

Figure 1 shows a time development of the double Maxwell–Jüttner distribution through the collisional relaxation. The double-peaked distribution function approximates to the equilibrium state whose velocity shift and temperature are $\mathbf{v}_{0}/2$ and $2\theta$, respectively. Figures 2(a) and (b) indicate the errors of conservation laws for the structure-preserving and intuitive discretizations, respectively. The structure-preserving scheme with the collision kernel of Eq. (17) strictly preserves the conservation laws of mass, momentum and energy only with round-off errors. It seems that the mass, momentum of $z$-direction and energy slightly accumulate the round-off errors. We reported that Padé-type (or implicit) filters used in computational fluid dynamics can accumulate round-off errors, but it can be suppressed by changing the order of arithmetic operations Shiroto et al. (2017). Therefore, the accumulation of round-off errors can be resolved by an optimization of the code although it is out of scope of this article. From the other perspective, the conservation property of the proposed scheme is quite fine as well as it is affected by the order of arithmetic operations. In contrast, the intuitive discretization with the collision kernel of Eq. (7) clearly violates the energy conservation. Note that the mass and momentum are strictly conserved with the intuitive scheme. This is because the relativistic LFP equation is discretized in the divergence-form, and the intuitive collision kernel also satisfies the symmetry of $\mathsf{U}(\mathbf{u},\mathbf{u^{\prime}})=\mathsf{U}(\mathbf{u^{\prime}},\mathbf{u})$.

Finally, we address some remaining issues about the positivity and H-theorem. Figure 3 shows the H-theorem is strictly preserved and the linear relaxation rate is initially reproduced well in the current situation. Here we introduced the entropy as follows to avoid negative values due to round-off errors:

[TABLE]

The linear relaxation rate is obtained as follows by assuming the initial Maxwellian:

[TABLE]

Reference Berezin et al. (1987), for example, suggests that there is an upper-bound of $\Delta t$ which ensures the positivity and H-theorem. In this experiment, the collision time is resolved quite well so that the requirement is satisfied. However, such a small temporal interval is not suit for practical simulations. We did not give a discussion about the temporal discretization in this article since the conservation laws do not depend on the temporal structure of the relativistic LFP equation. A development of conservative and entropic scheme for the relativistic LFP equation will be performed in the separate paper.

VI CONCLUSIONS

A feasibility of conservative scheme for the relativistic Landau–Fokker–Planck equation has been demonstrated. The proposed scheme has a unique way of calculating the collision kernel specialized for linear finite-difference operators. The verification via thermal-equilibration problem manifests the conservation of mass, momentum and energy only with round-off errors. Although there are still some problems of computational cost, positivity, H-theorem and boundary conditions, our strategy gives a piece of puzzle for practical simulation of the collisional relaxation in the relativistic regime.

Acknowledgements.

T.S. would like to appreciate valuable information from Dr. William T. Taitano and Dr. Luis Chacón (Los Alamos National Laboratory). This work was supported by KAKENHI Grant Numbers JP15K21767 and JP18H05851. Numerical experiments were carried out on NEC SX-ACE, Cybermedia Center, Osaka University.