Particle Physics Playground

Four forces in a box

Every particle is a disc with a mass $m$ and a charge $q$, moving in the plane of the screen. Positions are in pixels and times in seconds, and the charges and masses are in arbitrary units, so the numbers on the sliders are the numbers in the equations. Two particles feel Coulomb's electric force [1] and Newton's gravity [2], combined into one pair force along the line between them,

where $\mathbf r_{ij}$ points from $i$ to $j$ and $\mathbf F_{ij}$ acts on $j$. Like charges repel, unlike charges attract, and gravity always attracts. The softening length $\varepsilon$, 12 pixels by default, is Plummer's trick [4], standard in gravitational simulations [5]. It caps the force when two particles pass through each other, which would otherwise be infinite, and leaves it unchanged at large distances. A uniform magnetic field $B_z$ perpendicular to the screen adds the Lorentz force [3], which bends the path of every moving charge without touching its speed,

With the field pointing out of the screen, positive charges circle clockwise and negative ones counterclockwise. A fourth, optional force is a linear drag that slowly brakes everything. Random particles get masses between 0.5 and 2.5, charges of size 6 to 12 with either sign, and velocity components up to 30 pixels per second. A disc's radius grows with its mass as $3\sqrt{m}+2$ pixels.

A step that keeps the energy

The simulation advances in fixed steps of $\Delta t = 1/480$ s of simulated time, whatever the frame rate of the screen. The Speed slider changes how many steps run per frame, never the step itself, so the motion is the same on every device. Each step is a half kick from the electric and gravitational forces, a rotation of every velocity by the angle $qB_z\Delta t/m$, a second half kick, and a drift of the positions. This kick, rotate, kick sequence is the Boris scheme [6], the standard way to move charges through magnetic fields in plasma simulations, and it owes its success to preserving phase-space volume [7]. Here the rotation uses the exact angle, so the magnetic field changes no speed at all.

Without a magnetic field the scheme reduces to the leapfrog or Störmer–Verlet method. It is symplectic: it does not conserve the energy exactly, but its error stays bounded instead of growing with time [8]. The adaptive Runge–Kutta integrator of the three-body page is more accurate per step but not symplectic, so its small energy error can creep in one direction. The measured numbers show it. With walls, no losses and elastic or no collisions, the total energy of the Coulomb Chaos, Mini Collider and Plasma in B presets stays within three parts per million of its starting value over 30 simulated seconds. All numbers quoted on this page were measured with the engine itself, in a box of 800 by 600 pixels. Each step costs a pass over all pairs, about 31,000 of them for 250 particles.

Leapfrog hides one subtlety. It stores velocities half a step out of phase with the positions, so the kinetic and potential energies of the readout must be evaluated at the same instant, with the velocity $\mathbf w = \mathbf v + \mathbf a\,\Delta t/2$. Walls and collisions act on the same $\mathbf w$. They reverse its normal component and never move a particle, so a lossless bounce leaves the total energy exactly as it was.

Bounces, merges and edges

In an elastic collision two overlapping discs that approach each other exchange an impulse along the line of their centres, which reverses their relative velocity in that direction and conserves both momentum and kinetic energy. In a merging collision they become one disc at their centre of mass, with their combined mass, charge and momentum. The kinetic energy of their relative motion is lost, which is how clumps grow in the Gravity N-Body preset. Its 70 neutral bodies merge to 35 within 10 seconds, to 24 by 30 seconds, and to 7 after about 100 seconds. Switching the lossless option off brings in the drag and makes the walls return only 90% of the normal speed, and the energy then drains away.

The edges need not be walls. On a torus a particle that leaves on the right comes back on the left, and on a Möbius strip it also comes back upside down, mirrored top to bottom. The forces do not reach across the seam, though, so the energy jumps a little whenever a particle wraps around.

Charges in a magnetic field

A charge in a uniform magnetic field circles at the cyclotron frequency $\omega_c = |q|B/m$, on a circle of radius $v/\omega_c$ [9]. In the Plasma in B preset, with $B_z = 1.8$, a typical particle goes round once every 0.58 seconds on a circle with a radius of only two or three pixels, smaller than the particle itself. The particles are tied to their place, and that changes the meaning of the electric force. A charge held in a strong field does not move along the force but drifts at right angles to it, with the drift velocity $\mathbf v_d = \mathbf E\times\mathbf B/B^2$, the same for positive and negative charges [9].

Two particles show what this means. Released at rest 40 pixels apart without a field, a positive and a negative charge rush together. With $B_z = 1.8$ they keep their distance and glide side by side, at right angles to the line joining them, at 0.71 pixels per second, the value $E/B$ predicts. Two like charges, which would fly apart, instead circle each other at a steady 2 degrees per second. Each charge is carried by the field of the others, exactly as point vortices in a two-dimensional fluid carry one another. This analogy between a strongly magnetized plasma and a fluid of vortices is a classic of plasma physics [10].

The field can also alternate, $B_z(t) = B_0\sin 2\pi f t$, which reverses the sense of the circling twice per period. Nature would add something here. By Faraday's law a changing magnetic field induces an electric field, and that electric field does work. It is how the betatron accelerates electrons [11]. The playground leaves the induced field out, so its alternating field still changes no energy. Measured at 0.5, 2 and 5 Hz, the total energy stays within ten parts per million over 30 seconds.

Inside the code

The heart of the engine is the step itself. The accelerations come from one pass over all pairs, with each pair force computed once and applied to both particles with opposite signs.

// One fixed step of 1/480 s: half kick, magnetic turn, half kick, drift.
computeAccelerations(ps, s);              // Coulomb + gravity at the current positions
for (const p of ps) { p.vx += p.ax * half;  p.vy += p.ay * half; }
rotateInField(ps, fieldAt(s, t), dt);     // turn every velocity by q B dt / m
for (const p of ps) {
  p.vx = (p.vx + p.ax * half) * damp;     // damp = 1 unless drag is switched on
  p.vy = (p.vy + p.ay * half) * damp;
  p.x += p.vx * dt;
  p.y += p.vy * dt;
}

The magnetic part is a pure rotation of each velocity. Approximating it with a small force kick would let the speed creep up step after step, and a charge in a strong field would slowly spiral outward.

// The magnetic force only turns the velocity. Turning it by the exact
// angle keeps every speed, so the field never does work.
const ang = (p.q / p.m) * B * dt;
const c = Math.cos(ang), sn = Math.sin(ang);
const vx = p.vx;
p.vx = c * vx - sn * p.vy;
p.vy = sn * vx + c * p.vy;

Walls work on the synchronized velocity, so that a reflection conserves exactly the energy that the readout shows. Elastic collisions follow the same rule.

// Leapfrog keeps the velocity half a step out of phase with the position.
// The velocity at the same instant is w = v + a dt/2. A wall reverses w
// and leaves the position alone, so a bounce cannot change K + U.
const wx = p.vx + p.ax * half;
if ((p.x < r && wx < 0) || (p.x > width - r && wx > 0)) p.vx = -e * wx - p.ax * half;

Try it yourself

Press Clear, choose a positive charge under Charge of launched particles, and set the magnetic field to its maximum. Launch one particle with a long drag and watch it circle. Then press Clear again and launch a positive and a negative particle close together with very short drags, and see the pair drift off slowly side by side instead of colliding. Plain drags launch the charge you selected, and with a mouse a right drag launches a positive particle and a Shift or middle drag a neutral one. In the Mini Collider preset, with a weaker field, the circles are large enough to see, and fast particles curl more widely than slow ones. In Gravity N-Body, watch the bodies merge into a few heavy clumps and the total energy fall with every merger. Switch the edges to the Möbius strip and follow a particle around twice before it comes back the right way up. The readout of the change in total energy tells you at every moment whether the settings you chose conserve it.

References

1. [1] C.-A. Coulomb, Premier mémoire sur l’électricité et le magnétisme, Histoire de l’Académie Royale des Sciences, 569 (1785).
2. [2] I. Newton, Philosophiæ Naturalis Principia Mathematica (London, 1687).
3. [3] H. A. Lorentz, La théorie électromagnétique de Maxwell et son application aux corps mouvants, Arch. Néerl. Sci. Exactes Nat. 25, 363 (1892).
4. [4] H. C. Plummer, On the problem of distribution in globular star clusters, Mon. Not. R. Astron. Soc. 71, 460 (1911).
5. [5] S. J. Aarseth, Gravitational N-Body Simulations (Cambridge University Press, Cambridge, 2003).
6. [6] J. P. Boris, Relativistic plasma simulation: optimization of a hybrid code, in Proc. Fourth Conf. Numerical Simulation of Plasmas (Naval Research Laboratory, Washington, DC, 1970), p. 3.
7. [7] H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, and W. M. Tang, Why is Boris algorithm so good?, Phys. Plasmas 20, 084503 (2013).
8. [8] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration illustrated by the Störmer–Verlet method, Acta Numer. 12, 399 (2003).
9. [9] J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, New York, 1999).
10. [10] J. B. Taylor and B. McNamara, Plasma diffusion in two dimensions, Phys. Fluids 14, 1492 (1971).
11. [11] D. W. Kerst, The acceleration of electrons by magnetic induction, Phys. Rev. 60, 47 (1941).