Development of Coulomb Collision Operator for Multi-Ion Species Plasmas

Abstract

Understanding the interactions of multi-ion species is essential for analyzing particle transport in divertor plasmas. This study focuses on the development of a Coulomb collision operator utilizing the Nanbu collision algorithm for modeling Coulomb interactions among multi-ion plasmas composed of hydrogen isotope ions and electrons. Simulations initialized with a Maxwellian velocity distribution examine the energy relaxations and momentum exchanges among particle species. The simulation results confirm energy conservation and reveal energy relaxation patterns, with equilibrium timescales influenced by interspecies mass differences. Incorporating particles with unequal statistical weights improves computational efficiency, reducing simulation time while maintaining the accuracy of these relaxation dynamics. This approach accelerates the Coulomb interactions in multi-ion species plasmas, making it an effective tool in fusion plasma research.

1.  Introduction

Collisions play a crucial role in governing complex behaviors of the multi-ion species plasmas, especially in the divertor region of magnetic confinement fusion devices. The particle transport from the Scrape-Off Layer (SOL) to the divertor has created the multi-ion species plasmas, which consist of ions with a variety of masses and charges [1]. For the divertor, this region serves as a crucial part between the fusion plasma and the exhaust system, where ionized particles are removed or recycled [2]. In this region, transport of multi-ion species plasmas leads to complex interaction dynamics through collision processes. The collisions of the multi-ion species plasma are essential for enabling particles to exchange energy and momentum, as well as for controlling the flow of heat and particles to the divertor targets [3–5]. These processes help to support stable fusion core conditions, support plasma detachment from the divertor surface, and reduce the heat load on divertor components. The interpretation of collisions in multi-ion species plasmas becomes a key challenge due to the complex interactions related to mass, charge, and energy exchange among different species.

Due to the complex behaviors in multi-ion plasmas resulting from particle collisions, fully kinetic simulations are essential tools for accurately capturing physical interactions [6]. Unlike fluid or hybrid models, kinetic approaches resolve the phase-space dynamics of individual particles [7]. This allows for the accurate modeling of particle transport. Among kinetic methods, the Particle-In-Cell (PIC) simulation method [8] is a widely used tool for simulating collisionless and collisional plasmas. The PIC simulations track individual motions of particle species under self-consistent electromagnetic fields, and incorporate collisional models through the implementation of collision operators [9].

In the divertor of the fusion device, where the plasma density is relatively high and temperatures are lower than in the fusion core, particle collisions in the multi-ion species plasmas become dominant features. Collisions of the multi-ion species plasmas are primarily governed by Coulomb interactions, which are long-range forces arising from electrostatic interactions between charged particles [10]. These interactions typically involve frequent, cumulative small-angle scatterings. Accurately resolving these small scattering-angle results is a significant computational challenge, especially in large-scale simulations. To overcome this computational challenge, efficient numerical methods, such as binary collision algorithms [11, 12], have been developed to approximate the cumulative impact of small-angle collisions while significantly reducing the computational cost. One such binary collision algorithm is the Nanbu collision algorithm [13], which approximates the cumulative effects of multiple small scattering angles in the Coulomb collisions by grouping them into a single larger-angle event over a simulation time step. This collision algorithm significantly reduces computational cost while preserving essential features of collisional dynamics. The Nanbu collision algorithm can be further enhanced using a weighted particle method [14]. This allows for larger time steps and enables simulations of larger simulation systems or longer physical timescales. However, no existing work has specifically applied a Nanbu collision operator for multi-ion species plasma.

This study presents the development of a Coulomb collision operator based on the Nanbu collision algorithm for multi-ion plasmas. The collision operator accounts for Coulomb interactions between electrons and hydrogen isotope ions, such as protons and deuterium ions. The weighted particle method of the Nanbu collision algorithm is applied to particle collisions to observe energy relaxation features and reduce computational time. The energy relaxation features of the particle collisions are used to assess the accuracy of the simulations through the analysis of energy relaxation features from theoretical methods. Moreover, the energy relaxation features obtained using the weighted particle method are analyzed to assess the accuracy of the collisional system and the computational efficiency.

The rest of this paper is organized as follows: Sec. 1 introduces the study. Section 2 describes the methodology of the Coulomb collision operator, including the Nanbu collision algorithm, weighted particle method, and energy relaxation analysis. Section 3 presents the simulation parameters. Section 4 discusses the results, and Section 5 concludes the paper. An unnumbered Acknowledgment section follows.

2.  Methodology of Coulomb Collision Operator in Multi-Ion Species Plasmas

2.1  Method of Nanbu’s collision algorithm

Multi-ion species plasmas consist of ion species with different charges and masses that move freely in the numerical domain. The particle velocity distribution of each species follows a Maxwellian distribution, $f_s({\mathit{v}}_s)$, defined as:

  

$$f_s({\mathit{v}}_s)=n_s{\left(\frac{m_s}{2\pi k_{\mathrm{B}}{T}_s}\right)}^{3/2}\exp \left[-\frac{m_s{\mathit{v}}_s^2}{2k_{\mathrm{B}}{T}_s}\right],$$(1)

where $n_s$, $m_s$, ${\mathit{v}}_s$ and $T_s$ denote the density, the mass, the velocity, and the temperature of particle sth species, respectively. The Boltzmann constant is denoted as $k_{\mathrm{B}}$. The effects of Coulomb interactions perturb the particle distribution function and then introduce a Coulomb collision operator, $C_s(f_s(\mu ))$. The collision operator is expressed with a time derivative of the particle distribution function that: $C_s(f_s(\mu ))\equiv \partial f_s(\mu )/\partial s^{\alpha |\beta }$, where $f_s(\mu )$, $\theta$, $\varphi$, $s^{\alpha |\beta }$ are the particle distribution function in spherical unit, the polar angle, the azimuthal angle, and the normalized time in collisions of $\alpha -\beta$ species, respectively. The particle distribution function in the spherical unit is defined as:

  

$$f_s(\mu )\equiv {\int }_0^{\mathrm{\infty }}{\int }_0^{2\pi }{f}_s({\mathit{v}}_s){\mathit{v}}_s^{2}d\mathit{v}d\varphi ,$$(2)

where $\mu \equiv \cos \theta$. The normalized time, $s^{\alpha |\beta }$, is defined by $s^{\alpha |\beta }=\tau /{\tau }_{\text{coll}}^{\alpha |\beta }$, where $\tau$ and ${\tau }_{\text{coll}}^{\alpha |\beta }$ denote the physical time step size and the characteristic collision time, respectively. This dimensionless time parameter plays a crucial role in determining the evolution governed by the collision operator.

The Nanbu collision algorithm focuses on a single large-scattering angle, $\chi$, which represents the cumulative effect of multiple small-angle collisions, $\theta$. The primary goal of the collision algorithm is to achieve the accuracy of the collision process in plasma by accumulating small scattering angles. The Nanbu collision algorithm has slightly changed the probability function from the previous study [11] as a function of the scattering angle coupling with normalized time that: $f_s(\chi )=A\exp (A\cos \chi )/4\pi \sinh (A)$, where $A$ is the shape factor of the distribution function. The shape factor, $A$, is solved by the nonlinear equation: $\coth (A)=\exp (-s^{\alpha |\beta })+1/A$. The accumulative scattering angle, $\chi$, can be determined from the solution of $A$ that: $\cos \chi =(1/A)(\ln (e^{\mathit{-A}}+2U\sinh A))$, where $U$ is a random number generated from the Monte-Carlo method.

The normalized time, $s^{\alpha |\beta }$, according to the Nanbu collision algorithm is defined as:

  

$$s^{\alpha |\beta }=\frac{\ln {\mathrm{\Lambda }}^{\alpha |\beta }}{4\pi }{\left(\frac{q_{\alpha }{q}_{\beta }}{{ϵ}_0{\mu }_{\alpha \beta }}\right)}^2{n}_{\beta }{|\mathbf{g}|}^{\mathrm{-3}}\tau ,$$(3)

where $\ln {\mathrm{\Lambda }}^{\alpha |\beta }$, $q_{\alpha }$, $q_{\beta }$, ${ϵ}_0$, ${\mu }_{\alpha \beta }$, $n_{\beta }$, and $\left|\mathbf{g}\right|$ denote the Coulomb logarithm, the charge of the colliding species, the charge of the target species, the electrical permittivity in vacuum, the reduced mass of the collision species, the target density, and the magnitude of the relative velocity, respectively. The relative velocity is determined between the colliding species and the target species that:

  

$$\mathbf{g}\equiv {\mathbf{v}}_{\alpha }-{\mathbf{v}}_{\beta }.$$(4)

Finally, the post-collision velocities for colliding species, ${\mathbf{v}}_{\alpha }^{\prime }$, and target species, ${\mathbf{v}}_{\beta }^{\prime }$, are calculated using momentum and energy conservation in the laboratory frame. These post-collision velocities are given by:

  

$$\begin{array}{cc}\hfill {\mathbf{v}}_{\alpha }^{\prime }& ={\mathbf{v}}_{\alpha }-\frac{m_{\beta }}{m_{\mathit{\alpha \beta }}}[\mathbf{g}(1-\cos \chi )+\mathbf{h}\sin \chi ],\hfill \\ \hfill {\mathbf{v}}_{\beta }^{\prime }& ={\mathbf{v}}_{\beta }+\frac{m_{\alpha }}{m_{\mathit{\alpha \beta }}}[\mathbf{g}(1-\cos \chi )+\mathbf{h}\sin \chi ],\hfill \end{array}$$(5)

where $m_{\alpha }$, $m_{\beta }$, $m_{\alpha \beta }$, and $\mathbf{h}$ denote the colliding mass, the target mass, the total mass in $\alpha -\beta$ collisions, and the random perpendicular vector orthogonal to the relative velocity, respectively. The components of $\mathbf{h}$ vector are expressed with Cartesian coordinates as:

  

(6)

which $ϵ\equiv 2\pi U$, and

  

(7)

  

(8)

Updating particle motion with the post-collision velocities is efficient in maintaining the accuracy of plasma dynamics.

2.2  Weighted particle method

The Nanbu collision algorithm also allows one to use a larger time step to maintain the precision of Coulomb interactions [14]. In this part, the Nanbu collision algorithm has introduced a weighted particle factor, $W_s$. The weighted particle factor is defined as: $W_s\equiv n_s/N_s$, where the number of superparticles and the density of each species are denoted as $N_s$ and $n_s$, respectively. On the basis of the realistic behaviors of the divertor region, multi-ion species plasmas are assumed to exhibit quasi-neutrality. This means that the contribution amounts of particle species in plasma composition are not of the same weight at all. For the simplest case, the weighted particle factors for each species in each collision are equal, the number of densities and the number of superparticles satisfy $n_{\alpha }=n_{\beta }\equiv n$, and $N_{\alpha }=N_{\beta }\equiv N$, leading to $W_{\alpha }=W_{\beta }\equiv W$. Consequently, each pair of superparticles interacts only once, with equal weighting applied to both particles. Under these conditions, the physical time step size, $\tau$, in Eq. (3) for the simplest case is given by:

  

$$\tau =\frac{N\times W{\tau }_{\mathrm{coll}}^{\alpha |\beta }}{n}.$$(9)

In the case where $N_{\alpha }\ne N_{\beta }$, this follows that $W_{\alpha }\ne W_{\beta }$, the physical time step size, $\tau$, the Eq. (3) is given by:

  

$$\tau =\{\begin{array}{cc}\frac{\max (W_{\alpha },W_{\beta }){\tau }_{\mathrm{coll}}^{\alpha |\beta }}{W_{\beta }},\hfill & \text{when}{N}_{\alpha }\ge N_{\beta },\hfill \\ \frac{\max (W_{\alpha },W_{\beta }){\tau }_{\mathrm{coll}}^{\alpha |\beta }}{W_{\alpha }},\hfill & \text{when}{N}_{\beta }>N_{\alpha }.\hfill \end{array}$$(10)

2.3  Energy relaxation of Coulomb collision

Energy relaxation due to the Coulomb collisions is assumed to drive the equilibration of each collision in the multi-ion species plasmas. The multi-ion species plasmas in this study are initialized with isotropic Maxwellian distributions for each species, with unequal initial temperatures across species. This reflects the natural thermal distribution of the particle species in the plasma source. Collisions drive the multi-ion species plasmas toward thermal equilibrium through the exchange of momentum and energy between interacting particle species. As particle species collide, they exchange kinetic energy and momentum through velocity interactions. This presents a relaxation process that gradually equalizes temperatures across species. The analytical expression of the energy relaxation process in different species collision can be expressed as:

  

$$\frac{\Delta T}{dt}=\frac{T_{\beta }-T_{\alpha }}{{\tau }_{\text{relax}}},$$(11)

where $\Delta T$, $T_{\alpha }$, $T_{\beta }$, $t$, and ${\tau }_{\text{relax}}$ denote the temperature changes in collision process, the temperatures of particle species $\alpha$ and species $\beta$, the physical time in collision process, and the energy relaxation time, respectively. The relaxation time is defined as the time scale over which the incident species exchanges energy with the target species through Coulomb collisions, based on their combined thermal velocity [15]:

  

$${\tau }_{\text{relax,the}}^{\alpha |\beta }=\frac{3{\pi }^{3/2}{ϵ}_0^2{m}_{\mathrm{\alpha }}^2{\left(v_{\text{th},\alpha }^2+v_{\text{th},\beta }^2\right)}^{3/2}}{n_{\beta }{q}_{\alpha }^2{q}_{\beta }^2\ln {\mathrm{\Lambda }}_{\alpha \beta }(1+m_{\alpha }/m_{\beta })},$$(12)

where $v_{\mathrm{th},\alpha }$ and $v_{\mathrm{th},\beta }$ are the thermal velocities of particle species $\alpha$ and $\beta$, respectively.

3.  Numerical Parameters

This study considers the multi-ion species plasmas consisting of electrons (e⁻), protons (H⁺), and deuterium ions (D⁺). The particle density of each particle species, $n_s$, for the multi-ion species plasmas corresponds to a quasi-neutrality property of the divertor plasma. Therefore, the particle density of e⁻, H⁺, and D⁺ is $n_{{\mathrm{e}}^{\mathrm{-}}}$ =1 × 10¹⁸ m⁻³, $n_{{\mathrm{H}}^{\mathrm{+}}}$ = 0.5 × 10¹⁸ m⁻³, and $n_{{\mathrm{D}}^{\mathrm{+}}}$ = 0.5 × 10¹⁸ m⁻³, respectively. The simulations of the multi-ion species plasmas are divided into two groups with different temperatures at the initial time, as groups (I) and (II). For the multi-ion species plasma group (I), the temperatures of e⁻, H⁺, and D⁺ are $T_{{\mathrm{e}}^{\mathrm{-}},0}$ = 100, $T_{{\mathrm{H}}^{\mathrm{+}},0}$ = 50, and $T_{{\mathrm{D}}^{\mathrm{+}},0}$ = 100 eV, respectively. The multi-ion species plasma group (II) has temperatures of each species as $T_{{\mathrm{e}}^{\mathrm{-}},0}$ = 100, $T_{{\mathrm{H}}^{\mathrm{+}},0}$ = 100, and $T_{{\mathrm{D}}^{\mathrm{+}},0}$ = 50 eV. The Coulomb logarithm, $\ln {\mathrm{\Lambda }}^{\alpha |\beta }$, in each collision pair is derived from the initial temperatures and densities [16], and indicates in a range of 15–17. The physical time step size, $\tau$, is determined using Eqs. (9) and (10) based on the number of superparticles used in pairings. The characteristic collision time, ${\tau }_{\text{coll}}^{\alpha |\beta }$, for each collision is computed from the inverse of the collision frequency, $1/{\nu }_{\text{coll}}^{\alpha |\beta }$, as applied under the fusion plasma conditions [16]. The total number of iterations, $N_t$, performed in the present simulations is 20,000.

The simulation domain in this study is independent of spatial coordinates and electrostatic field interaction. The simulation domain employs free boundary conditions, with no specified source injection points or defined plasma source terminations. At the initial step, the multi-ion species plasmas are injected into the simulation domain with a constant particle flux and a velocity distribution governed by the Maxwellian distribution function. The multi-ion species plasmas thus move freely in the simulation domain and interact through the Coulomb interactions.

Since the computational cost associated with small scattering angles in Coulomb collisions [11] is exceedingly high, superparticles are used to represent groups of real particles. In this simulation study, one superparticle corresponds to 1 × 10⁴ real plasma particles, and the total number of superparticles per species is denoted by $N_s$. To accelerate the Coulomb collisions of the multi-ion species plasmas, a weighted particle factor, $W_s$, is used to enable a larger time step size in the three cases. In the first case, equal weights are assigned to all particle species by setting $N_{{\mathrm{e}}^{\mathrm{-}}}$ = 5 × 10⁷, $N_{{\mathrm{H}}^{\mathrm{+}}}$ = 2.5 × 10⁷, and $N_{{\mathrm{D}}^{\mathrm{+}}}$ = 2.5 × 10⁷ superparticles, respectively. This results in equally weighted particle factors as: $W_{{\mathrm{e}}^{\mathrm{-}}}=W_{{\mathrm{H}}^{\mathrm{+}}}=W_{{\mathrm{D}}^{\mathrm{+}}}$, where $W_{{\mathrm{e}}^{\mathrm{-}}}$, $W_{{\mathrm{H}}^{+}}$, and $W_{{\mathrm{D}}^{\mathrm{+}}}$ denote the weighted particle factors of e⁻, H⁺, and D⁺, respectively. In the second case, unequal weights are used with setting $W_{{\mathrm{H}}^{\mathrm{+}}}=1.67W_{{\mathrm{e}}^{\mathrm{-}}}$, and $W_{{\mathrm{D}}^{\mathrm{+}}}=0.72W_{{\mathrm{e}}^{\mathrm{-}}}$, when we set $N_{{\mathrm{H}}^{\mathrm{+}}}$ = 1.5 × 10⁷, $N_{{\mathrm{D}}^{\mathrm{+}}}$ = 3.5 × 10⁷, and $N_{{\mathrm{e}}^{\mathrm{-}}}$ = 5 × 10⁷ superparticles. The third case also has unequal weights, where $N_{{\mathrm{H}}^{\mathrm{+}}}$ = 3.5 × 10⁷, $N_{{\mathrm{D}}^{\mathrm{+}}}$ = 1.5 × 10⁷, and $N_{{\mathrm{e}}^{\mathrm{-}}}$ = 5 × 10⁷ superparticles. Therefore, $W_{{\mathrm{H}}^{\mathrm{+}}}=0.72W_{{\mathrm{e}}^{\mathrm{-}}}$, and $W_{{\mathrm{D}}^{\mathrm{+}}}=1.67W_{{\mathrm{e}}^{\mathrm{-}}}$. In all cases, the combination of the superparticle numbers and their corresponding weight factors ensures that the quasi-neutrality condition is satisfied, meaning that the total electron density equals the sum of the ion densities.

4.  Simulation Results

4.1  Typical features of energy relaxation from collisions in multi-ion species plasmas

The typical energy relaxation behaviors of the multi-ion species plasma groups (I) and (II) in this section are derived from simulations using equally weighted particles, where $N_{\mathrm{e}}^{-}$ = 5 × 10⁷, $N_{\mathrm{H}}^{+}$ = 2.5 × 10⁷, $N_{\mathrm{D}}^{+}$ = 2.5 × 10⁷ superparticles, and $W_{{\mathrm{e}}^{\mathrm{-}}}=W_{{\mathrm{H}}^{\mathrm{+}}}=W_{{\mathrm{D}}^{\mathrm{+}}}$. Each superparticle within a given species has the same and uniform particle weight. When collisions are performed, we follow a systematic pairing strategy. For electron-ion collisions in each multi-ion species plasma group, each superparticle from the ion species is randomly paired with some superparticles from electrons within the plasma group using the Monte Carlo method. This ensures that energy exchange is statistically accurate while maintaining computational efficiency. For ion–ion collisions between different species, each superparticle of one ion species is randomly paired once with a superparticle from the other ion species, ensuring symmetric treatment in cross-species interactions.

The typical energy relaxation features of the multi-ion species plasma groups (I) and (II) are shown as a function of physical time in simulations, $N_t\tau$, in Fig. 1. In the overall equilibration process of the simulations shown in Fig. 1, the two groups of the multi-ion species plasmas are initialized with a Maxwellian velocity distribution for each particle species. The Coulomb interactions drive the system toward equilibration through momentum and energy exchange, ultimately approaching an equilibrium temperature, $T_{\mathrm{\infty }}$, which is given by:

  

$$\begin{array}{cc}\hfill T_{\mathrm{\infty }}& =\frac{{\sum }_s{n}_s{T}_{s,0}}{{\sum }_s{n}_s}\hfill \\ & =\frac{n_{{\mathrm{e}}^{\mathrm{-}}}{T}_{{\mathrm{e}}^{\mathrm{-}},0}+n_{{\mathrm{H}}^{\mathrm{+}}}{T}_{{\mathrm{H}}^{\mathrm{+}},0}+n_{{\mathrm{D}}^{\mathrm{+}}}{T}_{{\mathrm{D}}^{\mathrm{+}},0}}{n_{{\mathrm{e}}^{\mathrm{-}}}+n_{{\mathrm{H}}^{\mathrm{+}}}+n_{{\mathrm{D}}^{\mathrm{+}}}},\hfill \end{array}$$(13)

which was calculated to be 87.5 eV.

Fig. 1.  (A) Energy relaxation features of collisions from the multi-ion species plasma group (I): e⁻ (100 eV) + H⁺ (50 eV) + D⁺ (100 eV), and (B) from the multi-ion species plasma group (II): e⁻ (100 eV) + H⁺ (100 eV) + D⁺ (50 eV). Zoom-in panels highlight the statistical uncertainties (error bars) that arise from Monte Carlo sampling of the particle-collision pairs.

To better understand the contribution of each collision, we conducted further analysis of the simulations by isolating individual collision pairs, as summarized in Table 1. Unlike Fig. 1, which illustrates only system-level equilibration, Table 1 reveals the distinct energy relaxation rates associated with each collision pair, showing how these rates are modified by the presence of the third particle species. We examine the energy relaxations of each collision for the plasma group (I), as shown in Table 1. Among the collisions in the plasma group (I), the fastest energy relaxation occurs in e⁻ (100 eV)–H⁺ (50 eV), which reaches equilibrium in approximately $\approx \mathrm{12,932}({\tau }_{\text{relax,the}}^{{\mathrm{e}}^{\mathrm{-}}|{\mathrm{H}}^{\mathrm{+}}})$ to balance momentum and energy transfer. While the e⁻ (100 eV)–D⁺ (100 eV) collision that reaches equilibrium in approximately $\approx \mathrm{16,389}({\tau }_{\text{relax,the}}^{{\mathrm{e}}^{\mathrm{-}}|{\mathrm{D}}^{\mathrm{+}}})$. The slowest energy relaxation of this plasma group is from the collision D⁺ (100 eV)–H⁺ (50 eV), which reaches equilibrium in approximately $\approx 670({\tau }_{\text{relax,the}}^{{\mathrm{D}}^{\mathrm{+}}|{\mathrm{H}}^{\mathrm{+}}})$.

Table 1. Pairwise energy relaxation times are measured from simulations, ${\tau }_{\text{relax,sim}}^{\alpha |\beta }$, are compared with the theoretical energy relaxation time, ${\tau }_{\text{relax,the}}^{\alpha |\beta }$. The theoretical relaxation times, ${\tau }_{\text{relax,the}}^{\alpha |\beta }$, are calculated from the Eq. (12). Each ${\tau }_{\text{relax,sim}}^{\alpha |\beta }$ value corresponds to a specific collision pair and is measured under conditions where unlike-particle and like-particle collisions are included simultaneously. As a result, each relaxation time reflects the influence of concurrent collisions with other particle species.

Collision Types	${\tau }_{\text{relax,the}}^{\alpha |\beta }$ (s)	${\tau }_{\text{relax,sim}}^{\alpha |\beta }$ (s)
Multi-Ion Species Plasmas: Group (I)
e− (100 eV)–H+ (50 eV)	2.621 × 10−5	0.338
e− (100 eV)–D+ (100 eV)	2.621 × 10−5	0.430
D+ (100 eV)–H+ (50 eV)	7.783 × 10−4	0.522
Multi-Ion Species Plasmas: Group (II)
e− (100 eV)–H+ (100 eV)	2.621 × 10−5	0.441
e− (100 eV)–D+ (50 eV)	2.621 × 10−5	0.492
H+ (100 eV)–D+ (50 eV)	7.673 × 10−4	0.551

The collisions of each particle pair for the multi-ion species plasma group (II) are shown also in Table 1. The fastest collision is from e⁻ (100 eV)–H⁺ (100 eV), which reaches equilibrium in approximately $\approx \mathrm{16,809}({\tau }_{\text{relax,the}}^{{\mathrm{e}}^{\mathrm{-}}|{\mathrm{H}}^{\mathrm{+}}})$. The collision e⁻ (100 eV)–D⁺ (50 eV) in this group reaches equilibrium in approximately $\approx \mathrm{18,774}({\tau }_{\mathrm{relax,the}}^{{\mathrm{e}}^{\mathrm{-}}|{\mathrm{D}}^{\mathrm{+}}})$. The slowest energy relaxation in this plasma group is from the collision H⁺ (100 eV)–D⁺ (50 eV) that reaches equilibrium in approximately $\approx 718({\tau }_{\text{relax,the}}^{{\mathrm{H}}^{\mathrm{+}}|{\mathrm{D}}^{\mathrm{+}}})$. Across both multi-ion species plasma groups, the electron–ion collisions involving lighter ions equilibrate more rapidly than those involving heavier ions, reflecting the mass dependence of ion species in Coulomb interactions. In contrast, ion–ion collisions between different species show energy relaxation patterns, indicating slower equilibration compared to electron–ion collisions, due to the larger masses involved.

After observing the energy relaxation behaviors in the multi-ion species plasma groups (I) and (II) in Table 1, we found that the relaxation times are longer than the theoretical collision times associated with each interaction. We then compared the simulations for each plasma group with the analytical solutions derived from Eq. (11). These comparisons are shown in Fig. 2(A) for the H⁺ (50 eV)–D⁺ (100 eV) collision in plasma group (I), while the same collision in the plasma group (II) is shown in Fig. 2(B). The comparisons between the simulations and the analytical solutions yield an approximate expression for the relaxation time in our simulations as:

  

$${\tau }_{\text{relax,sim}}\propto \left(\frac{m_{\alpha }}{m_{\beta }}\right)\left(1+\frac{T_{\beta }}{T_{\alpha }}\right)\times {\tau }_{\text{relax,the}}^{\alpha |\beta }.$$(14)

This refers to the contributions of the masses and temperatures of the colliding species which influence the collisional processes that drive the system toward equilibration. The statistical measurements show good agreement with the analytical solutions obtained from the simulated energy relaxation processes, with associated errors of less than 10%.

Fig. 2.  (A) Comparison of analytical and simulated temperature evolution for the D⁺ (100 eV)–H⁺ (50 eV) collision in multi-ion species plasma group (I). (B) Comparison of analytical and simulated temperature evolution for the H⁺ (100 eV)–D⁺ (50 eV) collision in multi-ion species plasma group (II). Analytical solutions are shown as dotted lines; simulation results are shown as solid lines.

4.2  Energy relaxation features with weighted particle method

For the cases using the unequally weighted particle method, these cases are applied to the two plasma groups. In each of the multi-ion species plasma groups, all superparticles of a lower-grouped species are randomly matched with some superparticles of a higher-group species within the plasma group. The weighted particle factors are calculated for each particle species and uniformly assigned to all superparticles within that species. In the first case of unequal particle weights, we used $N_{{\mathrm{e}}^{\mathrm{-}}}$ = 5 × 10⁷, $N_{{\mathrm{H}}^{\mathrm{+}}}$ = 1.5 × 10⁷, and $N_{{\mathrm{D}}^{\mathrm{+}}}$ = 3.5 × 10⁷ superparticles in the two plasma groups, resulting in $W_{{\mathrm{H}}^{\mathrm{+}}}=1.67W_{{\mathrm{e}}^{\mathrm{-}}}$ and $W_{{\mathrm{D}}^{\mathrm{+}}}=0.72W_{{\mathrm{e}}^{\mathrm{-}}}$. For the second case of unequal particle weights, we set $W_{{\mathrm{H}}^{\mathrm{+}}}=0.72W_{{\mathrm{e}}^{\mathrm{-}}}$ and $W_{{\mathrm{D}}^{\mathrm{+}}}=1.67W_{{\mathrm{e}}^{\mathrm{-}}}$, when $N_{{\mathrm{e}}^{\mathrm{-}}}$ = 5 × 10⁷ superparticles, $N_{{\mathrm{H}}^{\mathrm{+}}}$ = 3.5 × 10⁷ superparticles, and $N_{{\mathrm{D}}^{\mathrm{+}}}$ = 1.5 × 10⁷ superparticles in both plasma groups. Then the physical time step sizes, $\tau$, for these simulation cases are determined by the Eq. (10).

For each collision involving weighted particle cases in plasma group (I) are shown in Figs. 3(A)–(C), while those for plasma group (II) are presented in Figs. 3(D)–(F). The solid lines represent cases in which equal weight factors are applied to the colliding species, while the dashed lines correspond to cases involving unequal weight factors. These results confirm that the use of unequally weighted particle factors preserves the accuracy of energy relaxation dynamics, yielding simulation outcomes consistent with those obtained using equally weighted particles. Moreover, this approach yields a computational speedup of approximately 45% by reducing the number of required superparticles without sacrificing physical fidelity. Thus, the unequal weighting method offers a reliable and efficient means of accelerating simulations of Coulomb interactions in multi-ion species plasmas.

Fig. 3.  Temperature evolution of the multi-ion species plasma group (I) (left column figures) are presented using weighted particle cases into the three particle species in collisions: e⁻–H⁺ (A), e⁻–D⁺ (B), and D⁺−H⁺ (C). While temperature evolution patterns of the multi-ion species group (II) (right column figures) using weighted particle cases are presented with the three particle species in collisions: e⁻–H⁺ (D), e⁻–D⁺ (E), and H⁺–D⁺ (F). In each collision, equally weighted particle case is indicated with solid lines, and unequally weighted particle cases are indicated with dotted lines.

5.  Conclusion

In conclusion, this study developed Coulomb collision operators for plasmas with multiple ion species using the Nanbu collision algorithm. The numerical results demonstrate energy relaxation patterns influenced by particle mass differences, with theoretical comparisons validating the simulation accuracy. Using unequal particle weights in collision simulations improves the Coulomb interactions and reduces the computational cost. The use of weighted particle methods enables acceleration of multi-ion species collision simulations by allowing larger time steps without sacrificing accuracy. This development can be integrated into future particle-in-cell (PIC) simulations to investigate plasmas with multiple ion species in the divertor region.

Acknowledgment

This work is performed on “Plasma Simulator” (NEC SX-Aurora TSUBASA) of NIFS with the support and under the auspices of the NIFS Collaboration Research program (NIFS22KISS005).

References

• [1]  P.C. Stangeby, The Plasma Boundary of Magnetic Fusion Devices, 1st ed. (CRC Press, Boca Raton, Florida, 2000).
• [2]  R.K. Janev, Atomic and Molecular Processes in Fusion Edge Plasmas, 1st ed. (Springer, New York, 1995).
• [3]   G.  Janeschitz et al., J. Nucl. Fusion 42, 14 (2002).
• [4]   V.P.  Budaev, J. Phys. At. Nucl. 79, 1137 (2016).
• [5]   S.I.  Krasheninnikov et al., J. Plasma Phys. 23, 055602 (2016).
• [6]   D.  Tskhakaya, J. Contrib. Plasma Phys. 52, 490 (2012).
• [7]   R.  Procassini et al., J. Nucl. Fusion 30, 2329 (1990).
• [8]   D.  Tskhakaya et al., J. Contrib. Plasma Phys. 47, 563 (2007).
• [9]   D.  Tskhakaya, Computational Many-Particle Physics, edited by  H.  Fehske et al. (Springer Berlin Heidelberg, Berlin, 2008), 161–189.
• [10]   L.  Spitzer, Physics of Fully Ionized Gases (Interscience Monograph Series, Interscience Publishers, 1962).
• [11]   T.  Takizuka and  H.  Abe, J. Comput. Phys. 25, 205 (1977).
• [12]   R.H.  Miller and  M.R.  Combi, Geophys. Res. Lett. 21, 16 (1994).
• [13]   K.  Nanbu, J. Phys. Rev. E. 55, 4642 (1997).
• [14]   K.  Nanbu and  S.  Yonemura, J. Comput. Phys. 145, 639 (1998).
• [15]   F.L.  Hinton, Handbook of Plasma Physics 1, (North-Holland Publishing Company, Amsterdam, 1983), Chapter 1.5.
• [16]   A.S.  Richardson and U.S. Office of Naval Research, 2019 NRL Plasma Formulary (2019), pp. 1–72.