Alternative pairwise nuclear fusion optimization algorithm in arbitrarily weighted particle-in-cell simulations

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Alternative pairwise nuclear fusion optimization algorithm in arbitrarily weighted particle-in-cell simulations

Qian Dong，

Chao-Zhi Li，

Chen Xie，

Zhi-Dong Chen，

De-Bin Zou，

Wen Luo ，

Tong-Pu Yu

Nuclear Science and Techniques

Vol.37, No.6

Article number 103

Published in print Jun 2026

Available online 25 Mar 2026

DOI：10.1007/s41365-026-01952-z

CSTR：32136.14.NST.2026.06103

12500

An optimization algorithm for pairwise nuclear fusion with arbitrarily weighted macroparticles is proposed and demonstrated using particle-in-cell (PIC) simulations. A faster computation of the Coulomb collision operator is achieved compared with that obtained through the classical nuclear fusion algorithm. Thus, the proposed algorithm can be implemented in the PIC simulation code with marginal effort. The algorithm is benchmarked in scenarios with like-particles, D(d,n)³He, unlike-particles, D(t,n)⁴He, thermonuclear fusion, and beam-target fusion. Moreover, the CPU time could be reduced significantly without reducing the simulation accuracy.

Particle-in-cell simulationNuclear fusionCoulomb collisionBinary collision

Introduction

Particle-in-cell (PIC) simulations provide insights into the kinetic properties of particles (such as velocity and position) during plasma interactions by tracking charged particles under the influence of both self-consistent and external electromagnetic fields. This technique provides an intuitive understanding of the kinetic behavior of field-particle interactions and is a vital tool for evaluating various plasma physics problems [1-3], such as inertial confinement fusion (ICF) [4, 5], ion acceleration [6, 7], laser–nuclear interactions [8], and positron generation [9-12]. Although PIC simulations of plasma physical processes are widely used, two critical challenges need to be addressed in practical applications: 1) numerical noise, which can be suppressed significantly with various methods [13] (e.g., variable weights [14], adaptive mesh refinement [15]), and 2) the more critical challenge of the high computational cost, which results in long simulation times even on parallel computers [16]. Therefore, optimizing PIC codes is important and necessary to increase the capabilities of the corresponding simulations [17, 18].

In PIC simulations, the kinetic modeling of plasmas is usually dominated by collective effects. Most PIC codes omit short-range Coulomb collisions between particles. At high temperatures ( keV) and relatively low plasma densities ( m^-3), the collision effects are usually considered marginal, and a collision-free approximation is valid in PIC simulations [19]. However, in plasmas with lower temperatures or significantly higher densities, the effect of short-range interactions between particles on the evolution of the systems needs to be addressed. Reducing the computational effort required to account for Coulomb collisions in PIC simulations has become a primary concern in the research on topics including ICF and strong-field ionization [20].

Initially, the particle distribution in plasma was determined by directly solving the Fokker–Planck kinetic equation. The key to this method is the calculation of the Rosenbluth potential to assess the contribution of Coulomb collisions to the scattering coefficient [21]. However, coupling this method with PIC simulations was initially complex. In addition, this method requires numerous simulated particles to accurately calculate the Rosenbluth potential. In practical PIC simulations, most models are based on the “binary collision” proposed by Takizuka and Abe (TA77) [22]. This model replaces direct collisions by calculating the Coulomb cross-section between particles in the same cell and randomly selecting the appropriate scattering angle. In addition, this methodology enables the simultaneous consideration of multiple elastic and inelastic collision channels in a pairwise process. Lavell et al. [23] applied the corrected binary collision model originally developed for charged-particle collisions to all the scattering channels considered. They validated the accuracy and utility of their model in a wide range of plasma environments. However, in binary collision modeling, the computational efficiency (i.e., computational speed and accuracy) of high-dimensional Monte Carlo (MC) high-dimensional integration methods depends strongly on the sampling method used. A simple and efficient MC method for computing arbitrary ion velocity distributions was proposed by Xie et al. [24]. A method similar to Monte Carlo collision (MCC) to simulate the transition region between collisionless plasma and Coulomb collision plasma is presented in Ref. [25]. A lattice-based “collision field” concept handles the collision process between different types of particles. This has significant advantages for high-dimensional simulations. Manheimer et al. [26] developed a Coulomb collision model for electron–electron and electron–ion collisions. It was simplified to an isotropic scattering model. The N97 [27] method uses the same order N binary-pairing procedure presented in TA77. However, it uses a different distribution function for the scattering angle. A grid-based binary model for Coulomb collisions was described by Cohen et al. [28]. It defines the unique particle counterpart of a particle on the grid by randomly sampling the velocity distribution. The model omits the process of direct pairing of particles in the simulation cell, significantly reducing the CPU runtime.

However, all these methods assume that the macroparticles have identical weights so that the real particle density is linearly proportional to the number of macroparticles in a cell. In many simulations, variable weights are necessary to render the computational task feasible. Using the appropriate pairing statistics, Miller and Combi [29] applied the TA77 method to differentially weighted particles. The algorithm preserves energy and momentum and produces suitable relaxation time scales compared with theoretical predictions. Subsequently, Nanbu and Yonemura [30] considered a Coulomb collision algorithm for weighted particles based on the cumulative properties of Coulomb collisions in plasma. The effectiveness of the proposed algorithm was verified using several examples. Higginson et al. [31] detailed a method for applying fusion to a PIC simulation with different particle weights and validated the method using analytical benchmarks in thermonuclear and beam-target fusion reactions. Unlike the work of Higginson et al. [31], Wu et al. presented an algorithm for pairwise fusion applicable to macroparticles of arbitrary weights at relativistic energies [32]. It is compatible with the Coulomb scattering algorithm widely used by Nanbu [30] and Sentoku [33] and can be implemented in the PIC simulation code with Coulomb scattering with marginal effort. Recently, Angus et al. [34] described a method for adjusting particle velocities post-scatter to restore the exact conservation of momentum and energy. They also illustrated its efficacy with various test problems.

In the present work, we described an efficient optimization method for the nuclear fusion algorithm and its implementation within the existing PIC code Smilei [35]. The aim was to develop collision-operator schemes with maximum efficiency while maintaining simulation accuracy. Based on the classical nuclear fusion algorithm described in Ref. [22], our algorithm simplifies particle pairing and reduces the computational overhead by eliminating the particle address randomization step. Collision particles are selected randomly. This significantly reduces the CPU time and improves the computational efficiency of the collision processes.

Optimization of the nuclear fusion algorithm

2.1

Theoretical background

The two fusion-reactant macroparticles a and b undergo fusion and create two products A and B with an energy gain of Q. This is described by Eq. (1): (1)In the center-of-momentum (CM) frame, a macrocollision may be considered as being associated with two interacting macroparticles and the local density of their respective species. In general, the fusion probability P of this interaction in the CM frame can be expressed as (2)where n_max is the maximum density of two particles a and b. v_rel is the relative velocity between the particles. The fusion reaction cross-section associated with the kinetic energy of the incident particle (denoted as σab) is expressed as follows: (3)Here, E_cm and E_lab denote the kinetic energies of incident particles in the CM and laboratory frames, respectively. represents the Gamow energy given by . Here, is the fine-structure constant, Zi denotes the atomic number, ħ represents the reduced Planck constant, and c is the speed of light. m_r is the reduced mass, as defined by . The exponential in Eq. (3) denotes the Coulomb barrier penetration probability. S(E_lab) is known as the astrophysical Sfactor. It varies gradually with the kinetic energy for most nuclear reactions. According to the parametric fitting formulae of Chulick et al. [36], the values of for the reactions D(d,n)³He and D(t,n)⁴He are 31.40 and 34.37, respectively. The cross-section and Sfactor are parameterized as functions of the energy of the reactant particle.

Note that the transformation of the coordinate system during the implementation of nuclear fusion in the code is mandatory. When a fusion reaction occurs, the fusion products are generated in the CM frame, transformed into the laboratory frame, and finally transformed back into the simulation frame.

In the laboratory frame, the coordinate origin is fixed at a point in the laboratory similar to the target nucleus in a nuclear experiment. Particle a moves at a relative velocity of v_rel. Therefore, the velocity in the simulation frame should first be transformed into the relative velocity in the laboratory frame by the transformation .

The sum of the kinetic energy of the incident particle a and target nucleus b relative to the center of mass is typically referred to as the kinetic energy of their relative motion. It is denoted as E_cm. It represents the kinetic energy of the CM frame. As shown in Fig. 1, the distance between the incident particle a and target nucleus b is denoted by x, whereas that between the center of mass and target nucleus is denoted by x_c. According to the definition of the center of mass, (4)Here, represents the velocity of incident particle a (denoted by va). Similarly, represents the velocity of the center of mass (denoted by v_cm). This equation indicates that particle a moves toward particle b with velocity va, whereas the center of mass moves toward particle b with velocity v_cm. Therefore, the velocity of a relative to the center of mass (or equivalently, the velocity of a in the CM frame) is given by (5)The velocity of b relative to the center of mass (denoted by vb in the CM frame) is equal in magnitude to the velocity v_cm of the center of mass relative to b, but opposite in direction. Hence, the following relationship is obtained: (6)The sum of the kinetic energies of particles a and b in the CM frame can be calculated using Eq. (5) and Eq. (6). (7)Subsequently, the fusion products are generated in the CM frame with the conservation of energy and momentum. The energy released in a reaction (generally referred to as the Q value) is given by the difference between the remaining mass energies of the reactants and products. Therefore, the total kinetic energy in the CM frame after the reaction is Q+E′. Here, Q is the energy released, and E′ represents the total kinetic energy of the product particles. This energy is then distributed among product particles A and B in proportion to their masses. In the CM frame, the total momentum of the particles is zero both before and after the reaction. This implies that . Therefore, the kinetic energy of the product particles can be expressed as follows (considering particle A as an example): (8)where mA and mB are the masses of the product particles. The velocity of particle A in the CM frame is (9)To facilitate the calculation of the fusion product, the momentum space is transformed such that the velocity of the center of mass is in the z-direction. This transformation is expressed as follows: (10)Here, θ represents the polar angle of the outgoing particle in the z-CM frame, and φ represents the azimuthal angle of the outgoing particle in this frame. In the CM framework, particle emission is assumed to be isotropic. This implies that it occurs uniformly in all directions with respect to polar angle θ. Azimuthal angle φ is selected from a uniform distribution between 0 and 2π. These angles are then used to calculate the velocity of the first product in the z-CM frame, as follows: (11)From the conservation of momentum , the velocity of the fusion product B in the CM frame is (12)The particle velocity vector should then be transformed back into the CM frame using the inverse of the previous transformation. (13) (14) (15)After the previous coordinate transformations and velocity calculations, the particle velocity vector is rotated back to the simulation frame using an inverse transformation. Finally, the obtained fusion product is propelled and considered as a conventional macroscopic particle in a simulation frame. Although this calculation is performed for all binary pairs in each simulation cell, most of the computational time required by the binary collision operator is spent randomizing and pairing particles rather than computing the collision kinematics.

Fig. 1

(Color online) Relationship between center-of-mass velocity and incident particle velocity

2.2

Optimizing the algorithm

Most PIC codes (including Smilei [35]) with different collision models [27, 37, 38] were inspired by TA77 [22]. The time steps of these codes follow a similar procedure, as follows:

(a)

Particle groups. The particles are stored as a linked list within each simulation cell. This data structure has significant advantages for parallel computations.

(b)

Particle address index randomization. This ensures the randomness of the subsequent particle pairing.

(c)

Pairing similar particles. Pairs of particles of the same species are determined from the top of the addresses.

(d)

Pairing of unlike-particles. Pairs of particles from different species are selected from the top of the addresses.

Finally, new particles are generated by the fusion of the paired particles. The original linked list is updated, and the product particle information is stored. The primary steps of the algorithm described above are summarized in Fig. 2.

Fig. 2

(Color online) Implementation process of the classical binary collision (TA77) and optimization algorithms (NEW). (a) Grouping of particles in the simulation cell. (b) Randomly altering the index order of particles within the simulation cell. (c) Pairing of like-particles of TA77. (d) Pairing of unlike-particles of TA77. (e) Pairing of like-particles of the optimization algorithm. (f) Pairing of unlike-particles of the optimization algorithm

The relaxation time used in the Coulomb collision calculation is more accurate than the standard approximation. However, it is computationally expensive and requires substantial CPU time. This presents a significant challenge for large-scale PIC simulations. To improve the computational efficiency while maintaining the accuracy of the Coulomb collision model, we implemented targeted optimizations using the classical binary collision algorithm.

First, we eliminated the particle address index randomization. Macroparticles may collide multiple times during the simulation period. However, the collision pair selection is, on average, performed once per execution of the optimization algorithm. This would be a significant time investment if particle address index randomization is considered at each time step. The selection of colliding particle pairs is achieved by random numbers. This ensures the randomness of particle pairing and eliminates the need for a particle address randomization operation.

We then optimized the pairing of like-particles. As illustrated in Fig. 2(e), when particle pairs belong to the same type, N_num denotes the total number of particles in the simulation cell. Two random numbers Ra and Rb are generated subsequently (). The indices of the paired particles are determined from and . This step is performed in a loop until the pairing of N_num/2 particles is completed. If and , pa is incremented by one. Otherwise, it is decremented by one. It is necessary to recursively apply the above steps to the remaining simulation cells until all the cells are paired. The pseudocode is provided in Algorithm 1.

Pairing of like-particles

Input: Information regarding all macroparticles in the simulation cell; Number of particles in the cell (N_num)

Output: Information of N_num/2 particle collision pairs

1 repeat

2 Generate two random ;

// Determine the index of particle pairs

3 , , if pa = pb

then

// If particle pairs have an equal index within valid bounds

4 pa = pa + 1;

5 else

// If particle pairs have an equal boundary index

6 pa = pa - 1;

7 end

8 end

9 until Select N_num/2 Pairs of colliding particles;

10 return N_num/2 particle colliding pairs;

Third, we optimized the pairing of unlike-particles. For different types of particle pairings, the total number of particles with a lower density Na and that of particles with a higher density Nb should be determined first. This is similar to the optimization of the pairing of like-particles, where two random numbers are generated and the indices of the pairs are determined by and . As shown in Fig. 2(f), the above procedure is repeated until Na particle pairing is completed. This process is then executed recursively across all the simulation cells until all the cells are particle-paired. The pseudocode is given in Algorithm 2.

Pairing of unlike-particles

Input: Information regarding all macroparticles in the cell; Number of particles with lower density in the cell (Na); Number of particles with higher density in the cell (Nb)

Output: Information of Na particle collision pairs

1 repeat

2 Generate two random ;

// Determine different density particle index

3 ;

4 ;

5 until Select N_a pairs of colliding particles;

6 return Na particle colliding pairs;

The spatial grid size in this procedure is of the order of the Debye radius. Because it is unnecessary to consider the Coulomb collisions between particles spaced farther apart than the Debye radius, the interactions between the particles in different simulation cells can be omitted. This also enables parallel computation of the simulation procedures. A flowchart comparison between the classic binary collision algorithm (TA77) and optimization algorithm (NEW) is shown in Fig. 3.

Fig. 3

(Color online) Flow chart comparison between the classical binary collision algorithm (TA77) and optimization algorithm (NEW)

Benchmarks

To evaluate the efficiency of the optimization algorithm, we present three typical simulations: D(d,n)³He and D(t,n)⁴He thermonuclear fusion and D(d,n)³He beam-target fusion. The results demonstrate that the optimization algorithm computes the Coulomb collision process with maximum efficiency without compromising physical accuracy. These benchmarks were implemented using the PIC code Smilei [35]. All the tests were performed on a disabled field. Simulations can use cross-sectional data derived from experimental measurements or parametric fitting equations. In this case, the second approach was preferred to ensure a direct comparison of the simulation results with the theoretical predictions.

It is important to emphasize that the optimization algorithm operates on probability-driven random pairing, rather than on the complete traversal of all particles. Although the number of pairings is limited by lower-density species, the optimization algorithm leverages MC methods to achieve a robust sampling of high-density species. In nuclear fusion, the velocity distribution of ions follows the Maxwell–Boltzmann distribution, wherein energetic particles (although sparse) contribute significantly to the fusion reaction [31]. Random pairing cumulatively covers the high-energy region using multiple independent samplings (cycling each time step) rather than relying on a single-step traversal. Consequently, the statistical rise and fall of random pairing is attenuated by time-averaging effects in the long-duration simulations. This ensures the accuracy of the kinetic evolution.

3.1

D(d,n)³He thermonuclear fusion

In the first benchmark, we created a simulation box with dimensions Lx = 40πc/ω_r and Ly = 40πc/ω_r. We divided it into a 100 × 100 grid, with each cell containing an equal number of D ions. The number density of D ions was 10²⁰ cm^-3, and the particles were weighted unequally. Thermonuclear fusion occurs when ions attain high temperatures. Because the cross-section for fusion depends strongly on the relative velocity of the particles, most fusion events occur between ions in the tails of the distribution, with kinetic energies several times higher than the mean energy of the plasma. Simulations were performed at different temperatures to evaluate the fusion reactivity and neutron spectra.

The neutron energy spectra obtained at different temperatures are shown in Fig. 4. TA77 and the optimization algorithm yielded similar neutron energy spectra. The full width at half maximum (FWHM) of the neutron energy spectra increased as the initial ion temperature did. To validate the accuracy of the energy spectra generated by the optimization algorithm, the number of particles per cell (PPC) was fixed at 1000 for all the simulations. As illustrated in Fig. 5(a), by varying the initial temperature of the D ions, we compared the effective FWHM of the simulated energy spectra with the analytical fit of Ballabio [39] using relativistic kinematics. The results reveal a remarkable agreement. The ultimate goal of the optimization algorithm is to improve the computational efficiency without altering the underlying physics. Figure 5(b) presents the CPU runtime for the D(d,n)³He thermonuclear fusion simulation with an initial temperature of 10 keV. Compared with the simulation results of TA77, the CPU runtime reduced by 30.5%. Although the improvement in computational efficiency may not be universal, it is more evident when simulating a larger number of particles in complex collisions.

Fig. 4

(Color online) Neutron energy spectra for D(d,n)³He fusion in plasmas with ion temperatures of (a) 5 keV, (b) 10 keV, and (c) 20 keV. Here, the simulation results with TA77 are shown in red and those with the optimized algorithm are shown in blue

Fig. 5

(Color online) (a) Neutron energy full width at half maximum (FWHM) simulated by the optimized algorithm (NEW) is compared with the FWHM from the work of Ballabio [39]. (b) Total computation time for the D(d, n)³He thermonuclear fusion simulation in both cases

3.2

D(t,n)⁴He thermonuclear fusion

To test the optimization algorithm in the case of unlike-particles and illustrate the case of unequal weighting, a simulation was conducted for D(t,n)⁴He thermonuclear fusion. The initial conditions were similar to those of the first benchmark, with a density of 10²⁰ cm^-3 for both ions and an initial temperature of 5 keV. The neutron energy spectra were benchmarked by varying the initial numbers of D and T ions in each cell. As shown in Fig. 6(a), when the initial number of D ions was less than that of the T ions (N_D = 100 ppc, N_T = 1000 ppc), the optimization algorithm matched the neutron energy spectrum from the TA77 simulation. In the other case, with the particle number reversed (N_D = 1000 ppc, N_T = 100 ppc), the results were as anticipated (see Fig. 6(b)). Furthermore, when equal weights were maintained and the initial density ratios of the D and T ions were varied, the neutron energy distribution was unaffected by the variations in initial particle density. It influenced only the peak of the energy spectrum. Because the nuclear reaction mechanism of the second benchmark is different from that of the first, we performed multiple sets of simulations. The evolution of the neutron energy FWHM with ion temperature is shown in Fig. 7(a). It again demonstrates the remarkable agreement between the simulation results of the optimization algorithm (NEW) and the theory (Ballabio). In addition, we compared the computational efficiency of D(t,n)⁴He thermonuclear fusion in both cases, as shown in Fig. 7(b)). This revealed a 25.5% reduction in the CPU runtime of the optimized algorithm compared with the TA77 results.

Fig. 6

(Color online) Neutron energy spectra for D(t,n)⁴He fusion in plasmas with ion temperatures of 5 keV. (a) Case where the initial number of particles is N_D = 100 ppc, N_T = 1000 ppc, respectively. (b) Case where the initial number of particles is N_D = 1000 ppc and N_T = 100 ppc. Here, the simulation results with TA77 are shown in red and those with the optimized algorithm are shown in blue

Fig. 7

(Color online) (a) Neutron energy FWHM simulated by the optimized algorithm (NEW) is compared with the FWHM from the work of Ballabio [39]. (b) Total computation time for the D(t,n)⁴He thermonuclear fusion simulation in both cases

3.3

D(d,n)³He beam-target fusion

Two distinct tests were conducted to investigate the feasibility of the D(d,n)³He beam-target fusion using an optimization algorithm. In the simulations, D(d,n)³He beam-target fusion occurred when a high-energy beam impacted a cold stationary D ions target. In this case, the kinetic energy of the beam determined the magnitude of the fusion cross-section.

One group of D ions was considered as a stationary background with a density of 10²⁰ cm^-3. The other group with an equal density was considered as a projectile. Two simulations with projectile velocities of 0.05 c and 0.1 c were conducted. The resulting neutron spectra are presented in Fig. 8. Here, Fig. 8(a) and Fig. 8(b) illustrate the simulation results with initial velocities of 0.05 c and 0.1 c, respectively. As shown, the optimization algorithm accurately reproduced the energy spectra of the neutrons. Because the optimization algorithm partially removed the particle traversal steps, we illustrate the relative error in the total kinetics during the simulation in Fig. 9(a). For the optimization algorithm, the peak of the relative error stabilized at ±0.02%. This is an insignificant energy deviation considering the effect of statistical fluctuations during the simulation. In addition, as shown in Fig. 9(b), the optimization algorithm significantly reduced the CPU runtime (by 32.5%).

Fig. 8

(Color online) (a) Neutron energy spectra for D(d,n)³He beam-target fusion when the initial projectile velocity is 0.05 c. (b) Neutron energy spectra for D(d,n)³He beam-target fusion when the initial projectile velocity is 0.1 c. Here, the simulation results with TA77 are shown in red and those with the optimized algorithm are shown in blue

Fig. 9

(Color online) (a) Relative error in the total kinetic energy owing to the use of the optimization algorithm during the simulation. (b) Total computation time for the D(d,n)³He beam-target fusion simulation in both cases

Conclusion

In this study, we focused on the computationally inefficient problem of multiparticle complex collisions in PIC simulations and optimized the classical PIC binary Coulomb collision algorithm. The particle-pairing process was simplified by removing the particle address randomization step. Additionally, collision particle pairs were selected using random numbers. This eliminated the need to traverse all the particle address information within the simulation cell and significantly reduced the CPU runtime. The optimized algorithm was benchmarked in like-particles D(d,n)³He, unlike-particles D(t,n)⁴He, thermonuclear fusion, and beam-target fusion. The results show that the CPU runtime can be generally reduced by approximately 30% while maintaining physical accuracy. The minimal errors in energy and momentum enable the optimized algorithm to minimize its effect on the underlying physics.

For plasmas with a significantly high density ( m^-3) or significantly low temperature ( keV), the collision frequency may increase significantly. This, in turn, would amplify the statistical error in the optimization algorithm (NEW). Under such conditions, the particle correlation effects (e.g., collective oscillations or strong coupling) may exceed the assumptions of the algorithm for independent collision events. We plan to further validate the robustness of the optimization algorithm in future studies. The generalization for parallel computations is apparent. This is because the optimization algorithm is performed cell-by-cell, and particles are “sorted” during the entire simulation. Therefore, the performance of the optimization algorithm is likely to improve substantially as the number of particles and complexity of the simulation system increase. More importantly, the optimization algorithm enables the evaluation of problems that would otherwise not be accessible owing to computational resource limitations. Thus, it provides technical support for subsequent large-scale particle collision simulations, including ICF and laser-driven neutron sources.

References

E.P. Alves, W.B. Mori, F. Fiuza et al.,

Numerical heating in particle-in-cell simulations with Monte Carlo binary collisions

. Phys. Rev. E 103, 013306 (2021). https://doi.org/10.1103/PhysRevE.103.013306

Alternative pairwise nuclear fusion optimization algorithm in arbitrarily weighted particle-in-cell simulations

Introduction

Optimization of the nuclear fusion algorithm

Theoretical background

Optimizing the algorithm

Benchmarks

D(d,n)3He thermonuclear fusion

D(t,n)4He thermonuclear fusion

D(d,n)3He beam-target fusion

Conclusion

D(d,n)³He thermonuclear fusion

D(t,n)⁴He thermonuclear fusion

D(d,n)³He beam-target fusion