NBODYCL: A Program to Generate Regularized Classical Trajectories for the n-Body Coulomb Problem

1. Introduction

The dynamics of charged particles interacting through the long-range Coulomb potential play a central role in many areas of physics, ranging from atomic and molecular collisions to plasma dynamics and astrophysical processes. However, the accurate numerical treatment of such systems remains a persistent challenge due to the singular nature of the Coulomb interaction at short inter-particle separations. In particular, direct integration of the classical equations of motion for few- and many-body Coulomb systems often leads to numerical instabilities or loss of energy conservation during close encounters of two or more particles. Overcoming these singularities is essential for constructing reliable trajectory-based simulations, specifically involving particles such as electrons and positrons.

Existing computational frameworks for charged-particle dynamics typically handle the Coulomb interaction either by introducing soft-core potentials, adaptive time steps, or cutoff distances to avoid singularities during close encounters. While such techniques can stabilize direct numerical integration, they do so at the cost of distorting the true Coulomb dynamics near short inter-particle separations. The true nature of the Coulomb potential at short distance is essential to accurately capture the dynamics in cases where, for example, an electron performs several encounters with an oppositely charged ion due to a laser field. It has been shown, that in laser-assisted collisions, the exact inclusion of the Coulomb potential between the electron and a proton leads to interesting phenomena [1,2,3,4,5] which are not observed if the Coulomb potential is softened. Some exact treatments of the Coulomb singularity include the treatment of frustrated ionization for strongly driven molecules in ref [6] and later analysis of double and triple ionization of three electron-atoms [7] and molecules [8], driven by infrared laser pulses.

In the present work, we introduce a general program to compute the classical trajectories in the $n-$charged particle system by globally regularizing the Coulomb singularities. It must be noted that the regularizing singular collisions in the $n-$body problem is not new, rather it has been extensively studied and successfully applied in the problem of gravitational $n-$body problem. See, for example, the classic texts by Aarseth [9] and Szebehely [10], and the references therein. More practical approach to computing regularized trajectories in the $n-$body gravitational problem is due to the works of Heggie [11] and Mikkola [12]. Therefore, we closely follow their derivations to arrive at the formula for the Coulomb system. In addition, we also introduce an oscillating electric field to the charged-particle system to address a larger class of problems in collisions physics.

The organization of the present paper is as follows. In Section 2 we derive our main equations, and in Section 3 we describe the main program nbodycl. For comparisons, we also provide a standard generalized integrator program called nbody-direct for the $n-$body Coulomb system without regularization. In Section 4 we compare the results from nbodycl and nbody-direct applied to example problems with number of charged particles $n=2,3,4$ and 5. Finally in Section 5 we discuss the additional challenges in the nbodycl program and its potential applications.

2. Formulation

For an n charged-particle system interacting with a time-dependent electric field, the system’s Hamiltonian $H^{\prime }$ in the dipole approximation can be written as,

$$\begin{array}{c}\hfill H^{\prime }=\sum _i^n{\displaystyle \frac{{{\mathbf{p}}^{\prime }}_i^2}{2m_i}}+\sum _i^{n-1}\sum _{j=i+1}^n{\displaystyle \frac{e_i{e}_j}{|{{\mathbf{q}}^{\prime }}_i-{{\mathbf{q}}^{\prime }}_j|}}-\sum _i^n{e}_i{{\mathbf{q}}^{\prime }}_i·\mathbf{F}\left(t\right),\end{array}$$

(1)

where the canonical variables ${{\mathbf{q}}^{\prime }}_i$ and ${{\mathbf{p}}^{\prime }}_i$ are the Cartesian coordinate and its conjugate momentum of the $i^{\mathrm{th}}$ particle. The mass and the charge of the particle are given by $m_i$ and $e_i$, respectively. The form of the time-dependent electric field vector $\mathbf{F}\left(t\right)$ can be arbitrary. Note that we will be using atomic units meaning that for an electron (positron) $m_i=1\left(1\right)$ and $e_i=-1(+1)$. The solution paths of Eq. (1) are given by the integration of Hamilton’s equations,

$$\begin{array}{cc}\hfill {{\dot{\mathbf{q}}}^{\prime }}_i& ={\displaystyle \frac{\partial H^{\prime }}{\partial {{\mathbf{p}}^{\prime }}_i}}={\displaystyle \frac{{{\mathbf{p}}^{\prime }}_i}{m_i}},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {{\dot{\mathbf{p}}}^{\prime }}_i& =-{\displaystyle \frac{\partial H^{\prime }}{\partial {{\mathbf{q}}^{\prime }}_i}}=\sum _{j\ne i}^n{e}_i{e}_j{\displaystyle \frac{({{\mathbf{q}}^{\prime }}_i-{{\mathbf{q}}^{\prime }}_j)}{|{{\mathbf{q}}^{\prime }}_i-{{\mathbf{q}}^{\prime }}_j{|}^3}}+e_i\mathbf{F}\left(t\right),\hfill \end{array}$$

(3)

where the dot represent the time derivative with respect to t. The direct numerical integration of the Hamilton’s equations (Newton’s equations) in Eqs. (2) and (3) introduce severe numerical instabilities when two or more particles approach infinitesimally close to each other, i.e., when $|{{\mathbf{q}}^{\prime }}_i-{{\mathbf{q}}^{\prime }}_j|\to 0$. These instabilities are commonly known as Coulomb singularities where the kinetic and potential energies of the system diverge. Therefore, we look for a global regularization for the Hamiltonian in Eq. (1) without Coulomb singularities.

For the regularization, we introduce the new coordinate variable ${\mathbf{q}}_{ij}$ and its conjugate momentum ${\mathbf{p}}_{ij}$ for each pair of $ij-$particles. For an n distinguishable particle system, there are $n(n-1)/2$ distinct pairs of particles. We can count the particle pairs with index k as the following code snippet shows. We also note that the indices $ij$ can be directly mapped to the index k via $k=k(i,j)=(i-1)n-i(i+1)/2+j$. The new coordinate variables are defined as follows,

$$\begin{array}{c}\hfill {\mathbf{q}}_k:={\mathbf{q}}_{ij}={{\mathbf{q}}^{\prime }}_i-{{\mathbf{q}}^{\prime }}_j,\end{array}$$

(4)

$$\begin{array}{c}\hfill \langle \mathbf{q}\rangle ={\displaystyle \frac{1}{M}}\sum _i^n{m}_i{{\mathbf{q}}^{\prime }}_i,\end{array}$$

(5)

where $M={\sum }_i^n{m}_i$ and $1\le i<j\le n$. The inverse transformation for the coordinate variable can be obtained from Eqs. (4) and (5) to give,

$$\begin{array}{c}\hfill {{\mathbf{q}}^{\prime }}_i=\langle \mathbf{q}\rangle +{\displaystyle \frac{1}{M}}\sum _{j=i+1}^n{m}_j{\mathbf{q}}_{ij}-{\displaystyle \frac{1}{M}}\sum _{j=1}^{i-1}{m}_j{\mathbf{q}}_{ji}.\end{array}$$

(6)

The new conjugate momenta are given by,

$$\begin{array}{c}\hfill {\mathbf{p}}_k:={\mathbf{p}}_{ij}={\displaystyle \frac{1}{n}}\left({{\mathbf{p}}^{\prime }}_i-{{\mathbf{p}}^{\prime }}_j-{\displaystyle \frac{m_i-m_j}{M}}\langle \mathbf{p}\rangle \right),\end{array}$$

(7)

$$\begin{array}{c}\hfill \langle \mathbf{p}\rangle =\sum _i^n{{\mathbf{p}}^{\prime }}_i.\end{array}$$

(8)

Note that if we wish to stay in the center-of-mass frame, we can take $\langle \mathbf{q}\rangle =0$ and $\langle \mathbf{p}\rangle =0$. As before, the inverse transformation for the conjugate momentum can be obtained from Eqs. (7) and () to give,

$$\begin{array}{c}\hfill {{\mathbf{p}}^{\prime }}_i=\sum _{j=i+1}^n{\mathbf{p}}_{ij}-\sum _{j=1}^{i-1}{\mathbf{p}}_{ji}+{\displaystyle \frac{m_i}{M}}\langle \mathbf{p}\rangle .\end{array}$$

(9)

The inverse transformation in Eq. (9) can be recast into the following compact form,

$$\begin{array}{c}\hfill {{\mathbf{p}}^{\prime }}_i=\sum _k^K{a}_{ik}{\mathbf{p}}_k+{\displaystyle \frac{m_i}{M}}\langle \mathbf{p}\rangle ,\end{array}$$

(10)

where $K=n(n-1)/2$, $a_{ik}=1$ and $a_{jk}=-1$ when $k=k(i,j)$, otherwise $a_{ik}=0$. In other words (recall that index k refers to the $(i,j)$ pair of particles with $i<j$), if the particle i is the first index of the $(i,j)$ pair, then $a_{ik}=1$, and if the particle j is the second index of the $(i,j)$ pair, then $a_{jk}=-1$, and finally, if the particle i is not included in the pair given by k, then $a_{ik}=0$. Following code snippet shows how the $a_{ik}$ terms are calculated. With the new sets of conjugate variables $[{\mathbf{q}}_k,\langle \mathbf{q}\rangle ]$ and $[{\mathbf{p}}_k,\langle \mathbf{p}\rangle ]$, and plugging Eqs. (6) and (10) in Eq. (1), the transformed Hamiltonian H can be written as,

$$\begin{array}{c}\hfill H=\sum _{k,k^{\prime }}^K{T}_{k{k}^{\prime }}{\mathbf{p}}_k·{\mathbf{p}}_{k^{\prime }}+{\displaystyle \frac{1}{2M}}{\langle \mathbf{p}\rangle }^2+\sum _k^K{\displaystyle \frac{e_k}{|{\mathbf{q}}_k|}}-\sum _k^K{\beta }_k{\mathbf{q}}_k·\mathbf{F}\left(t\right)-\sum _i^N{e}_i\langle \mathbf{q}\rangle ·\mathbf{F}\left(t\right),\end{array}$$

(11)

where we have defined the new symbols,

$$\begin{array}{c}\hfill T_{k{k}^{\prime }}={\displaystyle \frac{1}{2}}\sum _i^N{\displaystyle \frac{a_{ik}{a}_{i{k}^{\prime }}}{m_i}},e_k=e_i{e}_j,{\beta }_k=(e_i{m}_j-e_j{m}_i)/M\end{array}$$

(12)

The Coulomb singularities still exist in the transformed Hamiltonian H in Eq. (11) as the relative distances vectors ${\mathbf{q}}_k$ are still in the denominator. To remove the ${\mathbf{q}}_k$ from the denominator, we apply the transformation of variables proposed by Kustaanheimo and Stiefel (KS transformation) [9,13,14] and then later do a time transformation [11,12]. In the KS transformation, first the $K+1$ coordinate variables (all ${\mathbf{q}}_k$ plus $\langle \mathbf{q}\rangle$) are converted to new 4-vector variables ${\mathbf{Q}}_k=(Q_{k,1},Q_{k,2},Q_{k,3},Q_{k,4})$. Recall that the vector ${\mathbf{q}}_k$ is a 3-vector in coordinate space with cartesian components along $x,y$ and z directions. The relation between the ${\mathbf{q}}_k$ and ${\mathbf{Q}}_k$ is given below.

$$\begin{array}{cc}\hfill q_{k,1}& =Q_{k,1}^1-Q_{k,2}^2-Q_{k,3}^2+Q_{k,4}^2\hfill \\ \hfill q_{k,2}& =2(Q_{k,1}{Q}_{k,2}-Q_{k,3}{Q}_{k,4})\hfill \\ \hfill q_{k,3}& =2(Q_{k,1}{Q}_{k,3}+Q_{k,2}{Q}_{k,4})\hfill \\ \hfill q_{k,4}& =0\hfill \end{array}$$

(13)

The transformation in Eq. (13) can be cast into a matrix form ${\mathbf{q}}_k={\widehat{Q}}_k{\mathbf{Q}}_k$ where the square-matrix ${\widehat{Q}}_k$ is given by,

$$\begin{array}{c}\hfill {\widehat{Q}}_k=\left(\begin{array}{cccc}{Q}_{k,1}& -Q_{k,2}& -Q_{k,3}& Q_{k,4}\\ Q_{k,2}& Q_{k,1}& -Q_{k,4}& -Q_{k,3}\\ Q_{k,3}& Q_{k,4}& Q_{k,1}& Q_{k,2}\\ Q_{k,4}& -Q_{k,3}& Q_{k,2}& -Q_{k,1}\end{array}\right)\end{array}$$

(14)

Given ${\mathbf{q}}_k$, the construction of ${\mathbf{Q}}_k$ can be achieved as follows by defining $r=q_k$ and $u=\sqrt{(r+|q_{k,1}\left|\right)/2}$,

$$\begin{array}{c}\hfill {\mathbf{Q}}_k=\left(\begin{array}{c}u\\ q_{k,2}/u\\ q_{k,3}/u\\ 0\end{array}\right),\mathrm{if}{q}_{k,1}\ge 0\end{array}$$

(15)

$$\begin{array}{c}\hfill {\mathbf{Q}}_k=\left(\begin{array}{c}{q}_{k,2}/u\\ u\\ 0\\ q_{k,3}/u\end{array}\right),\mathrm{if}{q}_{k,1}<0\end{array}$$

(16)

Now, according to the KS-transformation, the relation between the momentum variable ${\mathbf{P}}_k$ conjugate to ${\mathbf{Q}}_k$ with the old momentum ${\mathbf{p}}_k$ is given by,

$$\begin{array}{cc}\hfill {\mathbf{p}}_k& ={\widehat{Q}}_k{\mathbf{P}}_k/\left(2Q_k^2\right),\hfill \\ \hfill {\mathbf{P}}_k& =2{\widehat{Q}}_k^T{\mathbf{p}}_k,\hfill \end{array}$$

(17)

where ${\widehat{Q}}_k^T$ is the transpose of the matrix ${\widehat{Q}}_k$ given in Eq. (14). With the new KS variables, the Hamiltonian H in Eq. (11) can be written as,

$$\begin{array}{cc}\hfill H(\mathbf{Q},\mathbf{P})& =\sum _{k,k^{\prime }}{\displaystyle \frac{T_{k{k}^{\prime }}}{4Q_k^2{Q}_{k^{\prime }}^2}}{\mathbf{P}}_k^T{\widehat{Q}}_k^T{\widehat{Q}}_{k^{\prime }}{\mathbf{P}}_{k^{\prime }}+{\displaystyle \frac{1}{8M{\langle Q\rangle }^2}}{\langle \mathbf{P}\rangle }^T{\langle \widehat{Q}\rangle }^T\langle \widehat{Q}\rangle \langle \mathbf{P}\rangle \hfill \\ \hfill & +\sum _k^K{\displaystyle \frac{e_k}{Q_k^2}}-\sum _k^K{\beta }_k{\mathbf{Q}}_k^T{\widehat{Q}}_k^T\mathbf{F}\left(t\right)-\sum _i^N{e}_i{\langle \mathbf{Q}\rangle }^T{\langle \widehat{Q}\rangle }^T\mathbf{F}\left(t\right).\hfill \end{array}$$

(18)

Now, based on the Hamiltonian $H(\mathbf{Q},\mathbf{P})$, we introduce the extended Hamiltonian K where we treat the time t as a new dynamical variable. Let us denote the canonical momentum to the time variable as $p_t$, then the extended Hamiltonian reads,

$$\begin{array}{c}\hfill K(\mathbf{Q},\mathbf{P},t,p_t)=H+p_t.\end{array}$$

(19)

With the above definition, it can be shown that $p_t\left(t\right)=-H\left(t\right)=-E\left(t\right)$ where $E\left(t\right)$ is the total energy along the solution path. As the final step in regularization, we employ a time transformation of the form,

$$\begin{array}{c}\hfill dt=g(\mathbf{Q},\mathbf{P})ds,\end{array}$$

(20)

where s is the pseudotime tracking the evolution of the paths. The time-transformation function $g(\mathbf{Q},\mathbf{P})$ is a matter of choice. The shapes of the trajectories in the configurations space should not depend on the choice for $g(\mathbf{Q},\mathbf{P})$. For example, it can be the product of all relative separations of the particles ${\Pi }_k^K{Q}_k^2$ [11], which eliminates all the $Q^2$ in Hamiltonian Eq. (18). In our calculations, we take $g(\mathbf{Q},\mathbf{P})=1/L$ where L is the Lagrangian of the system [12,15]. The system’s Lagrangian L behaves as $1/q_k$ when the separation between the $k^{th}$ particle pair approaches zero. Therefore, the function $g=1/L$ behaves as $q_k$ in this limit. After the time-transformation step, the final regularized Hamiltonian $\Gamma$ of the system is constructed as below,

$$\begin{array}{c}\hfill \Gamma (\mathbf{Q},\mathbf{P},s)=g(\mathbf{Q},\mathbf{P})K(\mathbf{Q},\mathbf{P},t,p_t;s).\end{array}$$

(21)

The regularized equations of motion can now be derived from Eq. (21). For this let us write, $L=\mathcal{T}-\mathcal{U}$, $H=\mathcal{T}+\mathcal{U}$, and $\Gamma =K/L$, where $\mathcal{T}$ and $\mathcal{U}$ are the kinetic and potential energy. Then, the Hamiltonian equations for $\Gamma$ can be written as,

$$\begin{array}{cc}\hfill Q_k^{\prime }={\displaystyle \frac{d{Q}_k}{ds}}& ={\displaystyle \frac{\partial \Gamma }{\partial P_k}}=g{\displaystyle \frac{\partial K}{\partial P_k}}+K{\displaystyle \frac{\partial g}{\partial P_k}}\hfill \\ \hfill & =g{\displaystyle \frac{\partial K}{\partial \mathcal{T}}}{\displaystyle \frac{\partial \mathcal{T}}{\partial P_k}}+K{\displaystyle \frac{\partial g}{\partial \mathcal{T}}}{\displaystyle \frac{\partial \mathcal{T}}{\partial P_k}}\hfill \\ \hfill & =g(1-\Gamma ){\displaystyle \frac{\partial \mathcal{T}}{\partial P_k}},\hfill \end{array}$$

(22)

and,

$$\begin{array}{cc}\hfill P_k^{\prime }={\displaystyle \frac{d{P}_k}{ds}}& =-{\displaystyle \frac{\partial \Gamma }{\partial Q_k}}=-g{\displaystyle \frac{\partial K}{\partial Q_k}}-K{\displaystyle \frac{\partial g}{\partial Q_k}}\hfill \\ \hfill & =-g\left[{\displaystyle \frac{\partial K}{\partial \mathcal{T}}}{\displaystyle \frac{\partial \mathcal{T}}{\partial Q_k}}+{\displaystyle \frac{\partial K}{\partial \mathcal{U}}}{\displaystyle \frac{\partial \mathcal{U}}{\partial Q_k}}\right]-K\left[{\displaystyle \frac{\partial g}{\partial \mathcal{T}}}{\displaystyle \frac{\partial \mathcal{T}}{\partial Q_k}}+{\displaystyle \frac{\partial g}{\partial \mathcal{U}}}{\displaystyle \frac{\partial \mathcal{U}}{\partial Q_k}}\right]\hfill \\ \hfill & =-g(1-\Gamma ){\displaystyle \frac{\partial \mathcal{T}}{\partial Q_k}}-g(1+\Gamma ){\displaystyle \frac{\partial \mathcal{U}}{\partial Q_k}}.\hfill \end{array}$$

(23)

The partial derivative of the kinetic and potential energies with respect to $Q_k$ and $P_k$ can be obtained from their forms found in Eq. (18) and they are also given in ref. [12] for the case of no external perturbation.

3. Program and Usage

3.1. Overview

We provide a Fortran 90 implementation nbodycl.f90, that integrates Coulomb $n-$body dynamics using a global KS regularization scheme as discussed in the previous section. The main subroutine takes in the particle details (masses, charges, and initial positions and momenta) and convert them to KS variables. It then advances all particle pairs (plus the center-of-mass) for a user-given time-step $dt$ before returning them to the main program in the original coordinates. The subroutine provided in the main program can be adapted to any general calculations where, for example, the user can change the initial positions and momenta randomly. An example driver program is provided within the main program for testing. All numerical integrations are performed using the Bulirsh-Stoer algorithm [16] with adaptive-step size control. The programs used in the present work are freely available at our Zenodo’s repository at the address; https://doi.org/10.5281/zenodo.17804482

3.2. Compilation and Execution

For compilation of the program, the intel Fortran 90 compiler is recommended. To compile, for example, use the command

ifort nbodycl.f90 -o nbodycl.exe

The inputs for the main subroutine can be generated inside the driver program or can be read from an external file. The example runs are prepared in this way and the structure of the input file is given below. nbodies is the number of charged particles in the system. The symbols m and q stand for the mass and charge of each particle.

nbodies

m1 q1

m2 q2

...

x1 y1 z1 px1 py1 pz1

x2 y2 z2 px2 py2 pz2

...

3.3. External Field (Optional)

The routine EFIELD(T,FLD) sets the dipole field and by default $\mathbf{F}\left(t\right)=0$. To drive the system with an oscillating electric field with a finite or long-pulse duration, the user can assign amplitudes, frequency,carrier envelope phase (CEP). The routine returns the Cartesian components of the field via FLD(1:3) array.

3.4. Limitations

The present version of the code can be used to obtain trajectories for upto $n\le 10$ bodies, but can be easily extended to treat $n>10$ by simple adjustments of the relevant dimensions of the arrays.

4. Results

In this section, we provide sample calculations for the case of $n=2,3,4$ and 5 body systems. To investigate the accuracy and advantage of the regularized method, we also compute the trajectories by directly integrating the Hamilton’s equations fom Eq. (1) which fail when two or more particles have close encounters. For this purpose, we also provide a direct integrator program nbody-direct.f90 which does not employ any regularization. In Figure 1 we show the example calculations for the cases where we have $n=2$ and 3 charged particles. They may as well represent the collisions between an electron and a proton, and a collision between a positron and a hydrogen atom, respectively. Similarly, in Figure 2 we show the examples cases with $n=4$ and 5 bodies. The error in energy is analyzed in both Figure 1 and Figure 2 and it is clear that our method significantly reduces such an error compared to direct integration.

5. Conclusions and Outlook

We have presented a global regularization scheme for the Coulomb N-body problem that combines pairwise Kustaanheimo–Stiefel (KS) coordinates transformation followed by a time transformation. The method removes close-encounter singularities globally without recourse to case-by-case treatments and advances the dynamics. At the end of the trajectory run, the program returns the laboratory coordinates for analysis. Our Fortran implementation (nbodycl) follows the derivation presented in the paper verbatim.

Presented results on the test cases with $n=2,3,4$ and 5 bodies demonstrate the smooth behavior of the trajectories through near-collisions and tight control of the total energy. The optional electric-field coupling enables laser-driven regularized Coulomb dynamics within the dipole-approximation. This makes the nbodycl program directly applicable to problems in AMO physics such as laser-assisted electron–ion scattering and positron/positronium interactions.