Authors: Oscar Eatwell, Jules Charlier, Jawad Ben Brahim, Rafael Salavarria Lorenzoni
This is an example screenshot of the collider at work:
Introduction
This is an interactive 3D visualization of electron-positron collisions, built entirely in Mathematica. Seeing as that future projects for colliders involve lepton colliders, we thought an electron-positron collider simulation would fit well, as it would be simulating what could be seen in the next decades. It is important to note that the main objective of this project is to visualise the collision, it is more of a "toy" simulation than a physically accurate one, as many parameters have been ignored for the sake of simplicity.
Code Explanation
Setup: Constants and Core Helpers
First, we need to set up our environment. This involves defining some fundamental physical and mathematical constants.
Clear[decayMemo];
decayMemo = <||>;
mMu = 0.105658; (** Muon mass in GeV/c^2 *)
tauMu = 2.196981*10^-6; (** Muon lifetime in seconds *)
cLight = 299792458.;
eps = 10^-6; (** A small number to prevent division by zero *)
Utility Functions
With the constants in place, we define a set of utility functions. These are the basic building blocks for everything that follows. They handle vector normalization, deterministic random direction generation (using BitXor on a seed), and a toy model for a three-body decay (DecayTriad) that conserves momentum in a simplified way. Momentum conservation is obviously a fundamental property we will be playing on in the 3D collision.
(** Normalizes a vector to unit length. Returns a default 'up' vector if the input vector has zero length. *)
unitVec[v_] := If[Norm[v] == 0, {0, 0, 1}, v/Norm[v]];
(** Generates a deterministic, pseudo-random 3D unit vector from an integer seed and an index. This ensures that the same event always looks the same. *)
fixedRandDir[seed_Integer, idx_Integer] :=
BlockRandom[SeedRandom[BitXor[seed, 1664525*idx + 1013904223]];
Module[{c, ph, s}, c = RandomReal[{-1, 1}];
ph = RandomReal[{0, 2 Pi}];
s = Sqrt[1 - c^2];
{s Cos[ph], s Sin[ph], c}]];
(** Converts spherical coordinate angles (theta, phi) to a Cartesian unit vector. *)
dirFromAngles[th_, ph_] := {Sin[th] Cos[ph], Sin[th] Sin[ph], Cos[th]};
(** Decomposes a parent particle's direction into three daughter directions, simulating momentum conservation in a toy model. *)
DecayTriad[uMu_, u1_, u2_, α_ : 0.35] :=
Module[{β = (1 - α)/2, ue, unu1, unu2},
ue = unitVec[u1];
unu1 = unitVec[u2];
(* Choose the third direction such that the weighted sum points along the parent direction uMu. *)
unu2 = unitVec[uMu - α ue - β unu1];
{ue, unu1, unu2}];
(** Safely calculates the end point of a particle track, ensuring it has a minimum length. *)
safeEnd[start_, dir_, L_, grow_] :=
start + (L*Max[grow, eps]) unitVec[dir];
(** A conditional graphics primitive. If 'cond' is True, it inserts 'expr' into the graphics list; otherwise, it inserts nothing. *)
gWhen[cond_, expr_] := If[TrueQ@cond, Sequence @@ {expr}, Sequence[]];
Graphics and Labeling
The next section of code is dedicated to functions that create and place labels for particles, handle special symbols (like $\nu$ for neutrino), and even includes a small visualization for the U(1) gauge symmetry of electromagnetism, represented by rotating "dials" as the beams approach. The goal is also to visualise the group theory behind the collision to make it more of an educative tool, since the mathematics behind the collision can sometimes seems abstract, our aim was to make it very visual.
(** ** Label and Glyph Helpers ** *)
(** Creates a styled 3D text label that uses TraditionalForm for mathematical expressions. *)
label3D[pos_, expr_, col_ : Black, sz_ : 14] :=
Inset[Style[expr, sz, Bold, col,
Background -> Directive[White, Opacity[.8]],
FormatType -> TraditionalForm], pos, Center, Scaled[.03]];
(** Generates the correct symbol for a neutrino (ν) or anti-neutrino (ν-bar). *)
nuGlyph[anti_ : False] := If[anti, Overscript["ν", "_"], "ν"];
gammaGlyph = "γ";
nuTauGlyph[anti_ : False] := If[anti, Overscript[Subscript["ν", "τ"], "_"], Subscript["ν", "τ"]];
chargeTag[sym_, sign_] := If[sign > 0, Row[{sym, "⁺"}], Row[{sym, "⁻"}]];
shortLabel[sym_, ch_] := Switch[sym, "gamma", "γ", "nu", "ν", "e", If[ch >= 0, "e⁺", "e⁻"], "mu", If[ch >= 0, "μ⁺", "μ⁻"], "pi", If[ch > 0, "π⁺", If[ch < 0, "π⁻", "π⁰"]], "rho", If[ch > 0, "ρ⁺", If[ch < 0, "ρ⁻", "ρ⁰"]], _, ToString[sym]];
(** Calculates a small, deterministic offset vector for a label to prevent it from overlapping with its particle track. *)
labelNudgeVec[u_, sid_, sym_, ch_] := Module[{a, v},
a = unitVec@fixedRandDir[sid, 4441 + Hash[{sym, ch}, "Adler32"]];
v = unitVec@perpComp[a, unitVec[u]];
If[Norm[v] < 10^-9, v = unitVec[{u[[2]], -u[[1]], 0}]]; v];
parComp[p_, u_] := (p.u) u;
perpComp[p_, u_] := p - parComp[p, u];
(** ** U(1) Gauge Symmetry Visualization ** *)
(** Creates an oscillating sphere to represent photon exchange during the beam approach. *)
drawPhotonExchange[p1_, p2_, tApproach_] :=
Module[{freq = 2 + 30*tApproach^2, phase = 2 Pi*tApproach, pos},
pos = p1 + (0.5 + 0.5*Sin[freq*phase])*(p2 - p1);
{Glow[Yellow], Yellow, Sphere[pos, 0.1]}];
U1Dial2D[phi_, col_ : Blue] := Graphics[{Thick, Circle[{0, 0}, 1], Arrow[{{0, 0}, 0.9 {Cos[phi], Sin[phi]}}], col, Disk[{0, 0}, .07]}, PlotRange -> 1.2, ImageSize -> 68];
U1Dial3D[pos_, phi_, col_ : Blue, scale_ : .06] := Inset[U1Dial2D[phi, col], pos, Center, Scaled[scale]];
drawPhoton[start_, dir_, L_, prog_] :=
Module[{u = unitVec[dir], end, pts, right},
end = start + L*Clip[prog, {0, 1}] u;
right = Normalize[Cross[u, {0.1, 0.8, 0.3}]]; (* A fixed axis for consistency *)
pts = Table[start + s (end - start) + (0.12 Sin[15 Pi s]) right, {s, 0, 1, 0.05}];
{Glow[Yellow], Yellow, Tube[pts, 0.03]}];
Particle Physics Engine: Decays and Channels
This is the core of the physics simulation. It contains tables for hadronic and leptonic decay modes (HadronicDecay, LeptonicDecay), functions to select an event type, and the branching ratios for complex decays like that of the tau lepton. It determines what exactly each particle will decay in based on probability and in which directions. This makes the collision simulation more physically meaningful, as it actually follows the experimental branching ratios of particles, creating decay chains.
(** A table of hadronic decay modes and their branching ratios. *)
HadronicDecay[sym_, ch_] :=
Switch[{sym, ch},
{"rho", 1}, {{1.0, {{"pi", 1}, {"pi", 0}}}},
{"rho", -1}, {{1.0, {{"pi", -1}, {"pi", 0}}}},
{"rho", 0}, {{1.0, {{"pi", 1}, {"pi", -1}}}},
{"pi", 0}, {{1.0, {{"gamma", 0}, {"gamma", 0}}}},
{"pi", 1}, {{0.999, {{"mu", 1}, {"nu", 0}}}, {0.001, {{"e", 1}, {"nu", 0}}}},
{"pi", -1}, {{0.999, {{"mu", -1}, {"nu", 0}}}, {0.001, {{"e", -1}, {"nu", 0}}}},
_, {}];
(** A table of leptonic decay modes. This handles the muon decays within the cascade framework. *)
LeptonicDecay[sym_, ch_] :=
Switch[{sym, ch},
{"mu", 1}, {{1.0, {{"e", 1}, {"nu", 0}, {"nu", 0}}}},
{"mu", -1}, {{1.0, {{"e", -1}, {"nu", 0}, {"nu", 0}}}},
_, {}];
weightedPick[list_] := With[{w = list[[All, 1]]}, RandomChoice[w -> list[[All, 2]]]];
(** ** Event Channel Selection ** *)
weightedChoice[assoc_Association] := With[{ks = Keys@assoc, ws = Values@assoc}, RandomChoice[ws -> ks]];
lep1Mix = <|"nunu" -> 0.200, "ee" -> 0.033, "mumu" -> 0.033, "tautau" -> 0.033, "gammagamma" -> 0.002|>;
highEMix = <|"nunu" -> 0.03, "ee" -> 0.03, "mumu" -> 0.03, "tautau" -> 0.02, "gammagamma" -> 0.01|>;
ChannelWeights[sqrtS_] := If[sqrtS < 170., lep1Mix, highEMix];
makeChannelState[sqrtS_] :=
Module[{ch = weightedChoice@ChannelWeights[sqrtS], seed = RandomInteger[{1, 10^9}]},
If[ch === "tautau",
<|"name" -> ch, "seed" -> seed, "tauModes" -> {weightedChoice@tauBR, weightedChoice@tauBR}|>,
<|"name" -> ch, "seed" -> seed|>]];
tauBR = <|"tau->e nu nu" -> 0.178, "tau->mu nu nu" -> 0.174, "tau->pi nu" -> 0.110, "tau->rho nu" -> 0.260, "tau->a1 nu(3pr)" -> 0.278|>;
The Decay Cascade Drawing Engine
To handle multi-stage decays (like $ \tau \rightarrow \rho + \nu_\tau \rightarrow \pi^+ +\pi^0 +\nu_\tau$, where the $\rho$ then decays to pions), we wrote a recursive drawing engine. decayCascadeDraw takes a particle and, if it's unstable, determines its decay products and then calls itself for each of them. This allows for arbitrarily deep decay chains. The animation is staged, so parent particles appear first, then their daughters, and so on. This is supposed to simulate the fascinating decay chains of particles, which can have many steps.
(** ** Deterministic Decay Choice and Geometry ** *)
a1ProngCharges[sign_, seed_] := Module[{key = {"a1_charges", seed, sign}, choice},
If[KeyExistsQ[decayMemo, key], Return[decayMemo[key]]];
choice = Block[{$RandomSeed = BitXor[seed, 884251]}, If[RandomReal[] < .5, sign*{1, 0, 0}, sign*{1, 1, -1}]];
decayMemo[key] = choice; choice];
childKey[sid_, idx_, depth_, {sym_, ch_}] := BitXor[sid, 73331 + 97*idx + 13*depth + 17*Hash[{sym, ch}, "Adler32"]];
decayChoice[sid_, idx_, depth_, sym_, ch_] := Module[{key = {sid, idx, depth, sym, ch}, modes, pick},
If[KeyExistsQ[decayMemo, key], Return[decayMemo[key]]];
modes = Join[HadronicDecay[sym, ch], LeptonicDecay[sym, ch]];
pick = If[modes === {}, {}, BlockRandom[SeedRandom[BitXor[sid, 73331 + idx + depth]]; weightedPick[modes]]];
If[ListQ[pick],
pick = MapIndexed[Append[#, First@#2] &, pick];
pick = SortBy[pick, {-Boole[#[[2]] != 0] &, #[[2]] &, childKey[sid, idx, depth, #[[;; 2]]] &, #[[3]] &}];
pick = pick[[All, ;; 2]];];
decayMemo[key] = pick];
clearDetRot[u_, sid_, key_, maxDeg_ : 22 Degree] := Module[{ax, th},
BlockRandom[SeedRandom[BitXor[sid, 1664525*key + 1013904223]];
ax = fixedRandDir[sid, key]; th = maxDeg*RandomReal[];
RotationMatrix[th, ax] . unitVec[u]]];
(** ** Core Drawing Functions ** *)
drawLeaf[start_, dir_, L_, sym_, ch_, labelSize_, showLabels_] := (* ... code from your input ... *);
decayCascadeDraw[start_, uParent_, L_, sym_, ch_, sid_, idxBase_, grow_, labelSize_, showLabels_, maxDepth_ : 4, depth_ : 0] := (* ... code from your input ... *);
drawTauStable[start_, sign_, seed_, mode_, eTrail_, grow_, labelSize_, showLabels_] := (* ... code from your input ... *);
(Note: The full code for drawLeaf, decayCascadeDraw, and drawTauStable is quite long, so I've omitted it here for brevity, but it's in the full notebook.)
The Main Scene
The makeScene function is the master controller. It takes all the parameters from the Manipulate controls (like beam energy, time, camera angle) and constructs the final Graphics3D object. It calculates the kinematics, determines the event channel, and then calls the appropriate drawing functions based on the event type.
(**
* This is the main function that constructs the entire 3D scene for a given event.
* It takes kinematics, viewing parameters, and event-specific controls as arguments.
*)
makeScene[
(** Kinematics *)
thN_?NumericQ, phN_?NumericQ, EN_?NumericQ, tN_?NumericQ,
(** Camera *)
viewDir_, viewUp_, camDist_?NumericQ,
(** Random Seeds for Muon Decay *)
eSeedMinus_, nuSeedMinus_, eSeedPlus_, nuSeedPlus_,
(** Animation and Display Controls *)
collisionEffectProgress_?NumericQ, showU1_, alphaU1_,
showDecays_, decayFrac_, eFrac_, eTrail_, showLabels_, labelSize_,
(** Event Selection *)
channelSel_ : "Auto"
] := (* ... full implementation from your input ... *)
The Interactive Manipulate Wrapper
Finally, the entire simulation is wrapped in a DynamicModule and Manipulate. This provides the user interface with sliders, buttons, and dropdown menus to control every aspect of the event. You can generate a new event, change the beam energy, play the animation, show or hide labels, and even choose a specific event channel to study, our favourite is the tautau channel, as it involves the most interesting decay chains.
SampleCosThetaQED[] := Module[{c, y}, While[True, c = RandomReal[{-1, 1}]; y = RandomReal[{0, 2}]; If[y <= (1 + c^2), Return[c]]]];
randDir[] := With[{c = RandomReal[{-1, 1}], ph = RandomReal[{0, 2 Pi}]}, Module[{s = Sqrt[1 - c^2]}, {s Cos[ph], s Sin[ph], c}]];
DynamicModule[{
(** State variables for the interactive session *)
costh = SampleCosThetaQED[], phi = RandomReal[{0, 2 Pi}],
eSeedMinus = randDir[], nuSeedMinus = randDir[],
eSeedPlus = randDir[], nuSeedPlus = randDir[],
channelState = <|"name" -> "mumu", "seed" -> 1|>,
lastChannelPick = "Auto", soundPlayed = False
},
Manipulate[
(* ... full Manipulate implementation from your input ... *)
]
]
The full notebook is attached below so you can play around with the simulation.
Made by Oscar Eatwell, Jules Charlier, Jawad Ben Brahim and Rafael Salavarria Lorenzoni.
Acknowledgements
Many thanks to Alex Upellini, Mario Veruete and Quentin Rible for organising and teaching the Wolfram Language Summer School 2025, without whom the project would have been impossible. Many thanks also to all the developers at Wolfram who continuously work to make a better language.
Attachments:
