$N$-body linear force law allowing analytic solutions

Abstract

We present a pair-wise force law in a system of $N$ particles that produces analytic solutions for arbitrary number of particles, masses, and initial conditions. Each pair of particles interacts via a force that is proportional to the product of their masses and their separation distance, with the force directed radially. We show that, despite the $N$-body interaction, each particle behaves as if it interacts only with the center of mass of the system. This effective two-body interaction behaves as Hooke’s Law with a common frequency for all particles, with the familiar analytic solutions for the trajectories. With these analytic solutions, it is possible to efficiently simulate a collection of these particles and incorporate other external forces. As an example, we simulate the particles within an adiabatically expanding container and calculate pressure and temperature in both the attractive and repulsive cases.

1 Introduction

The general $N$-body problem has long been a challenge in analytical mechanics. Newtonian gravity has been of particular interest, where the three-body problem admits analytical solutions only in special cases [1, 2, 3, 4, 5, 6, 7, 8, 9] and interactions between a sphere and a hoop generates intricate orbital motions [10]. Even with numerical solutions, the chaotic nature of $N$-body problems means that arbitrary systems will be unstable, presenting pedagogical challenges for students interested in simulating such a system.

This paper presents an investigation of a pairwise force law with analytic solutions in an $N$-body system. Intriguingly, this force law is shown to be equivalent to a central force acting on each particle, and each individual particle exhibits motion independent of its immediate neighbors.

Given their mathematical simplicity, these solutions can be conveniently integrated into undergraduate courses on advanced mechanics, within units on harmonic oscillators or gravitational interactions. Moreover, due to their computational efficiency, these systems can be accurately simulated with large $N$ values on personal computers, thus rendering them suitable for inclusion in computational projects or courses on numerical methods.

The absence of short range interactions, while convenient for quantum systems based on similar interactions [11, 12, 13, 14, 15], presents a lack of dynamism in the classical case. For this reason, further applications of this force law are most intriguing when paired with an external interaction. In this paper, we explore a thermodynamic application with the system in an expanding container, and suggest follow-up student projects.

The paper is organized as follows. In Section 2 the force law is introduced, and equivalence of an $N$-body interaction to a 2-body interaction with the center of mass is shown. In Section 3, we present thermodynamic results for attractive and repulsive systems in an expanding container. Conclusions are provided in Section 4. Supplemental materials include a number of short python programs. These include one for plotting the orbits of an N-particle system given a specified set of mass and initial condition values based on the analytic trajectories; a similar program that runs a numerical simulation of a similar system directly from the pair-wise forces. The collision algorithm is presented in Appendix A, while ideas for possible research projects for advanced undergraduate students involving these two systems are discussed in Appendix B.

2 Force Law and General Solution

Bertrand’s theorem shows that the only central forces that guarantee closed orbits are the inverse square and linear force laws [16]. Inspired by this, we introduce a toy model of a linear gravity-like force

$$\mathbf{F}=-Jm_{1}m_{2}\mathbf{r},$$ (1)

where $J$ is a constant with units of N/kg²m and can be positive (attractive force) or negative (repulsive). The force is that on mass $m_{2}$ due to mass $m_{1}$, and $\mathbf{r}$ is the vector from $m_{1}$ to $m_{2}$. The force is gravity-like in that the force between particles is proportional to the mass of both particles. It is spring-like in that it is linear. In the attractive case, the force is equivalent to a system of masses connected pair-wise by springs with spring constants

$$k_{ij}=Jm_{i}m_{j}$$ (2)

for each pair, as illustrated in Fig. 1.

Figure 1: Left: a visualization of the $N$-body force interaction. The magnitudes of the effective spring constants $k_{ij}=Jm_{i}m_{j}$ are represented by the thickness of the lines. Right: a visualization of this $N$-body force reduced to individual forces interacting with the center of mass (shown as an empty circle). The new effective spring constants are $k_{i}=JMm_{i}$, where $M$ is the sum of all masses.

2.1 Solutions for point particles

We solve the $N$-body problem by showing that it is equivalent to $N$ one-body problems. Consider a collection of $N$ particles with total mass $M$. The equation of motion of particle $i$ is

$$m_{i}\ddot{\mathbf{r}}_{i}=-J\sum_{k\neq i}m_{i}m_{k}(\mathbf{r}_{i}-\mathbf{r}_{k}).$$ (3)

The contribution from $k=i$ is zero, so it may be included for convenience. Including this term, the sum becomes

$$m_{i}\ddot{\mathbf{r}}_{i} = -Jm_{i}\sum_{k}m_{k}(\mathbf{r}_{i}-\mathbf{r}_{k})$$ (4)

$$= -Jm_{i}\left(\mathbf{r}_{i}\sum_{k}m_{k}-\sum_{k}m_{k}\mathbf{r}_{k}\right)$$ (5)

The sum in the second term is recognizable from the definition of the center of mass, and the equation of motion becomes

$$m_{i}\ddot{\mathbf{r}}_{i}=-JMm_{i}(\mathbf{r}_{i}-\mathbf{R}_{\mathrm{cm}}),$$ (6)

where $M=\sum_{k}m_{k}$ is the total mass of the system, and $\mathbf{R}_{\mathrm{cm}}$ is the position of the center of mass, which need not be stationary. Thus, the particle behaves as if it is governed by Hooke’s Law centered at $\mathbf{R}_{\mathrm{cm}}$, with effective spring constant $k_{i}=JMm_{i},$ as illustrated in Figure 1.

The solutions are the harmonic oscillator solutions

$$\mathbf{r}_{i}=\mathbf{A}\cos(\omega t)+\mathbf{B}\sin(\omega t)+\mathbf{R}_{\mathrm{cm}},$$ (7)

where $\omega=\sqrt{JM}$ for all particles, and the constant vectors $\mathbf{A},\,\mathbf{B}$ are determined as usual from initial conditions. In the absence of external forces, the center of mass of the system moves at constant velocity, and we can use coordinates with $\mathbf{R}_{\mathrm{cm}}=0$ so that the particles move in closed elliptical orbits with a common frequency, consistent with Bertrand’s Theorem [3, 17, 18, 19].

Thus, the $N$-body interaction has an analytical solution, with the motion of each particle effectively described as an interaction with the center of mass of the system, and the trajectories are described by analytical expressions. This differs from the general $N$-body problem, which cannot be solved analytically. Additionally, particles in this system remain bound, in contrast to the inverse-square law, which can eject particles from the system.

2.2 Particle interacting with continuous mass distribution

We now generalize the above results to a particle interacting with a continuous mass distribution. Consider placing the point mass m at the origin, and it interacts with a mass $M$ distributed in space within a volume $V$, with a position-dependent density $\rho(\mathbf{r})$. In that case, the force on $m$ due to a small volume element $dV$ within the extended body is

$$d\mathbf{F}=-Jm\rho(\mathbf{r})rdV.$$ (8)

The total force is found by integrating over the volume of the extended body

$$\mathbf{F}=-Jm\int\rho(\mathbf{r})rdV=-JmM\mathbf{R}_{\mathrm{cm}}.$$ (9)

Notice, due to the linear nature of the force law, there is no concern of singularities in allowing the particle of interest to be located within the larger mass distribution. From (9), it is clear that the shape, size, and orientation of the mass distribution have no effect on the motion of test mass $m$. The distributed mass does not need to be a rigid body, and it may rotate in any arbitrary manner, as long as it is rotating about its center of mass, which is how unconstrained bodies rotate [3, 20, 21].

2.3 Potential energy

The potential energy $U$, of the system is also found analytically. The potential energy of a system of particles attached via the properly chosen springs is

$$U=\frac{1}{2}\left[\frac{1}{2}J\sum_{i,k}m_{i}m_{k}(\mathbf{r}_{i}-\mathbf{r}_{k})\cdot(\mathbf{r}_{i}-\mathbf{r}_{k})\right],$$ (10)

where the additional factor of $\frac{1}{2}$ corrects for double counting. The case $i=k$ does not contribute to the sum, so it need not be removed. Expanding the product,

$$U=\frac{J}{4}\left[\sum_{i,k}m_{i}m_{k}(r_{i}^{2}+r_{k}^{2}-2\mathbf{r}_{i}\cdot\mathbf{r}_{k})\right],$$ (11)

and separating the terms yields

$$U=\frac{J}{4}\left[\sum_{i}m_{i}\sum_{k}m_{k}r_{k}^{2}+\sum_{k}m_{k}\sum_{i}m_{i}r_{i}^{2}-2\sum_{i}m_{i}\mathbf{r}_{i}\cdot\sum_{k}m_{k}\mathbf{r}_{k}\right].$$ (12)

The last term is the expression for the location of the center of mass, which is set to be zero. The summations over mass alone in the first two terms gives the total mass $M$, so the first two terms are identical. The resulting potential energy

$$U=\sum_{k}\frac{1}{2}JMm_{k}r_{k}^{2}$$ (13)

is identical to $k$ independent harmonic oscillators with effective spring constant $JMm_{k}$.

Alternatively, we calculate the work to move the $N$ particles from the origin to their positions $\mathbf{R}_{i}$ from the force law (1) directly, with the assumption that the center of mass remains at the origin. The work done in moving a single particle is

$$W_{k}=\int_{0}^{R_{k}}JMm_{k}r_{k}dr_{k},$$ (14)

as the work is non-zero only for the radial component of displacement. Performing the integration and summing over all particles,

$$\sum W_{k}=\sum_{k}\frac{1}{2}JMm_{k}r_{k}^{2},$$ (15)

which is equal to the potential energy calculation (13).

2.4 Repulsive Case

For a repulsive interaction, the above analysis is repeated with $J<0$. The solutions become exponentials

$$\mathbf{r}_{i}=\mathbf{C}e^{\Omega t}+\mathbf{D}e^{-\Omega t}+\mathbf{R}_{\mathrm{cm}}$$ (16)

where $\Omega\equiv\sqrt{-JM}$. The energy analysis is the same as in Section 2.2, with $J\rightarrow-J$. Unless $\mathbf{C}=0$, the particle escapes to infinity, with exponentially increasing velocity. The potential energy of a particle is always negative

$$U_{k}=-JMm_{k}r_{k}^{2},$$ (17)

which explains why the kinetic energy grows without bound as the particle escapes to infinity.

3 Applications

The existence of analytic solutions allows for efficient simulation of large numbers of particles, effectively modeling a gas of strongly interacting particles. In this section, we use analytic simulations to study thermodynamic properties of such a large collection.

Numerical techniques explicitly evolve the system forward in time given the current state of positions and velocities. It is pedagogically useful to show that this numerical technique produces results that match the analytic solutions of Section 2. Animating the trajectories using VPython is also instructive, as it effectively illustrates the “non-interaction” of nearby particles. This also demonstrates the benefits of the Euler-Cromer method, which conserves energy in oscillatory motion, while the Euler method does not. Numerical techniques are also needed if the system is subjected to external forces. Examples of these applications are described in Appendix B.

Any numerical technique using the original pair-wise force law (1) must calculate the distance between each particle pair, a problem that scales as $\mathcal{O}(N^{2})$. Additionally, errors in the differential equation solver depend on the time step, with the Euler-Cromer method having errors $\mathcal{O}(dt)$. Thus, simulating a large number of particles interacting for a long time is computationally intensive for a student project. This motivates the use of the analytic solutions when possible.

In the case of an isolated system, the analytic solutions (7) and (16) give the position and velocity of any particle at arbitrary time. In the case of a system inside a container, a continuous collision algorithm is used to handle any particles that collide with the container during a given time step. Any particle colliding with the wall affects the rest of the masses, because the change in velocity of that single particle affects the center of mass. The effects of the collisions on the velocity of the center of mass are accumulated while the check for collisions is made. The analytic solutions are updated to account for the change in center of mass velocity. The effects of a moving wall are also incorporated. This collision algorithm is elaborated in Appendix A.

3.1 Thermodynamics

Here the attractive and repulsive systems are examined in a 3D spherical container. A spherical container is necessary because a rectangular prism will not thermalize the motion of the particles, as the motion along each of the Cartesian coordinates is independent. A collection of 4,000 particles is created with spherically symmetric initial positions and velocities. The masses are chosen at random to be between 4 and 5 base units.

(a)

(b)

(c)

Figure 2: This figure shows relationships between thermodynamic quantities for the attractive case. Temperature (a) and pressure (b) both decrease as volume is increased. While pressure and volume have a clearly-defined relationship, it does not follow a power law.

In contrast to the ideal gas model, systems with long-range interactions do not have well-defined thermodynamic quantities that apply to the entire system [22, 23, 24]. Instead, we examine the pressure and temperature at the boundary of the container. The pressure is calculated by summing over the momentum of the particles that hit the wall in each time step,

$$P=\frac{\sum|\Delta p_{r}|}{4\pi R^{2}\Delta t},$$ (18)

where $p_{r}$ is the radial component of the particle’s momentum R, is the radius of the container, and $\Delta t$ is the computational time step. The temperature is determined by averaging the kinetic energy of the particles that hit the wall in a given time interval. The pressure and temperature are therefore only defined in terms of the interaction with the environment via the walls of the container. This is similar to the definition of thermodynamic parameters for bound gravitational systems [22, 23, 24]. In the attractive case, pressure and temperature are well defined only when the container is sufficiently small that the particles interact with the walls of the container.

We change the container radius quickly, allowing the particles to undergo free expansion into the increased volume. The system is then given time after each free expansion to come back to equilibrium. The thermalization time is particularly long for the repulsive case, as expanding the container changes the kinetic energy of the system notably. The values presented in the Figures 2 and 3 are the values at the end of the time segment used to equilibrate. All quantities are scaled by their values after an initial thermalization period. It is worthwhile to note that the only quantities being held fixed in these simulations is the number of particles and the total energy of the system.

(a)

(b)

(c)

Figure 3: Relationships between thermodynamic quantities for the repulsive case. (a) Because the kinetic energy increases without bound, temperature increases with volume. (b) The repulsive force increases as the particles get farther apart, so pressure increases as the container’s volume increases.

In examining Figures 2-3, it is evident that there is more scatter in the temperature values than in the pressure values. This is because the pressure depends on the radial component of the velocity of the particles that hit the wall, which will equilibrate due to collisions with the spherical container. The temperature also depends on the tangential component of the velocity, which is unaffected by contact with the wall.

For the attractive case, the relationships among pressure, temperature, and volume share some similarity with those for the ideal gas model (see Figure 2). However, log-log plots reveal the relationships between pressure and volume and pressure and temperature do not follow power laws.

For the repulsive case, an increase in the volume allows the potential energy to grow very large and negative, so that the kinetic energy of the particles impacting the wall of the container increases. As such, Figure 3 shows that the pressure and the temperature both increase with increasing volume, contrary to the behavior of attractive systems and ideal gases. The relationship between pressure and volume appears linear on visual inspection of Figure 3(b), but a log-log fit reveals that $P\sim R^{2.5}$ for this data.

4 Conclusions

We present a pair-wise $N$-body interaction that is shown to be equivalent to $N$ one-body interactions with the center of mass of the system. In both attractive and repulsive case, the analytical trajectories are given by harmonic oscillator solutions. This is in stark contrast to other long-range interactions such Newtonian gravity, where only special case solutions exist for $N>2$.

It is counter-intuitive that the motion of each particle is independent of the location and motion of all of the other particles of the system, except through the common center of mass location. The particles all interact with each other, and yet experience no scattering. This behavior was confirmed with a numerical simulation of the system.

For the attractive system, the relationships found among the pressure, temperature, and volume are qualitatively similar to what one would expect, pressure and temperature decreasing with increasing volume, and reaching zero for some critical size, as the bound particles have an upper limit on their distance from the center of mass. A more detailed study comparing the thermodynamics of this system for small containers with that of the ideal gas, or van der Waals gas would be interesting. For small volumes, the system is expected to behave as free particles in a gas.

The mathematical form of the solutions, indicating the independent classical motion of each of the masses, implies a separable Hamiltonian. It has been known for some time that, with identical particles, the attractive case leads to a separable Hamiltonian. [13] The results presented here show that such a solution for the quantum equivalent [11, 12, 13, 14, 15] is also possible for N particles with arbitrary masses. The repulsive case is not bound, so there is not an equivalent quantum mechanical model system.

Appendix A Collision handling algorithm for particles in a container

This algorithm uses event-driven collision detection and elastic collision handling to simulate interactions between particles and the walls of a container. This is a discrete collision detection algorithm [25, 26], as opposed to the more computationally expensive continuous collision detection algorithm [27].

The algorithm begins by using the analytical trajectories to determine which particles will collide with the container wall in the next time step. This is done by checking if the future position of each particle, $r(t+\Delta t)$, exceeds the container’s radius $R$. For each of these particles, the collision is treated as elastic. The new velocity of the particle is computed on the assumption that it occurs at position $r(t)$ and time $t$.

The change in momentum for each particle due to the collision is calculated, and the total change in momentum is summed for all particles. This summed momentum change is used to update the velocity of the center of mass ($V_{\text{cm}}$). To this point, no particle positions have been changed.

Using the updated center of mass velocity and the current positions and velocities of the particles, new parameters for each particle’s trajectories are calculated. The positions and velocities of all particles are then advanced to $t+\Delta t$, ensuring that no further collisions occur within this time step. The algorithm then repeats this process for subsequent time steps.

As an example, the data provided involves $N=$ 40,000 particles over 500,000 time steps. This simulation runs on a standard laptop in less than 40 seconds. For scenarios involving fixed walls, or when expansion into a larger container is permitted, energy conservation is confirmed within computational precision.

Appendix B Suggested student projects

The thermodynamic application of Section 3 suggests several other investigations for a student researcher. Including elastic hard-sphere collisions would affect individual trajectories, but the state variables of the system would not be affected. This is similar to the model of the ideal gas. This might lead to thermalization on a shorter time scale. The container wall can be moved more quickly, allowing work to be done on the system. It would be interesting to explore possible hysteresis effects of a container moving inward and then allowed to expand.

Additionally, these systems can be simulated in situations that include outside forces. For example, linear drag could be applied to one or more of the particles, allowing students to illustrate known results, such as the ultimate dissipation of all energy in a system experiencing linear drag, and the return of the affected particle to its initial position [28]. One could also simulate an $N$-body attractive system released in a uniform gravitational field and allowed to bounce on a hard surface, tracking the motion of the center of mass and the characteristic size of the system as it experiences these outside forces.

A selection of glowscript programs is included [29], which simulate interactions with an external environment. Due to external interactions, these simulations rely on traditional numerical integration, rather than the analytic form of the solutions. As anticipated, these simulation run times scale as $N^{2}$, in contrast with the simulations that make use of the analytic solution, which scale linearly with $N$.