Exploring Particle Physics Simulations in Python

Editors note:

At New Math Data we strive to regularly push the envelope, test our skills, and generally innovate, in addition to a deep commitment to our more traditional practices such as analytics and data warehousing. It is in that spirit that we present this article, an experiment in data modeling, where we ask the question “can we model both the Standard Model of particle physics and elements of String Theory in an idiomatic language?”

Introduction

Particle physics, the cornerstone of our understanding of the universe’s fundamental structure, is a field that unravels the intricate behaviors and interactions of the smallest known particles. From quarks and gluons to electrons and neutrinos, the myriad particles governed by the Standard Model exhibit complex dynamics that are pivotal to the fabric of our reality. Simulating these interactions not only allows us to predict particle behavior but also to visualize phenomena that are otherwise beyond the reach of direct observation.

In this blog post, we embark on a journey to construct a particle physics simulation using Python. Our goal is to model the forces at play, simulate particle interactions such as decay, collisions, and scattering, and visualize their trajectories. By leveraging the power of computational simulations, we can gain insights into particle dynamics that enhance our theoretical and practical understanding of the subatomic world.

For on-going updates, follow our GitHub repo.

We will begin by laying the groundwork with the implementation of particle and field classes, followed by detailing the methods to simulate interactions. Next, we will delve into the dynamics of running the simulation and updating particle states. Finally, we will visualize the trajectories of these particles to better understand their interactions and behaviors. Join us as we explore the fascinating world of particle physics through the lens of computational simulations.

Building the Foundation – Particle and Field Classes

Defining Particles

At the heart of our simulation lies the Particle class, a blueprint for the various particles we aim to model. Each particle is characterized by fundamental properties such as mass, charge, spin, and velocity. These attributes are not just numerical values; they determine how particles interact with each other and with fields.

The Particle class encapsulates these properties and provides methods to compute derived quantities like energy and momentum. The mass of a particle, measured in MeV/c², dictates its resistance to acceleration. Charge, expressed in elementary charge units (e), determines how the particle responds to electric and magnetic fields. Spin, a quantum mechanical property, influences the particle’s statistical behavior and interaction dynamics.

class Particle:
    def __init__(
        self,
        name,
        symbol,
        mass,
        charge,
        spin,
        string_mode,
        decay_constant=0.0,
        decays_to=None,
        velocity=None,
    ):
        self.name = name
        self.symbol = symbol
        self.mass = mass
        self.charge = charge
        self.spin = spin
        self.string_mode = string_mode
        self.decay_constant = decay_constant
        self.decays_to = decays_to if decays_to else []
        self.velocity = (
            velocity
            if velocity
            else (
                random.uniform(-0.5 * c, 0.5 * c),
                random.uniform(-0.5 * c, 0.5 * c),
            )
        )
        self.energy = self.calculate_energy()
        self.momentum = self.calculate_momentum()
        self.positions = [(0, 0)]  # Store positions for visualization

    def __repr__(self):
        return f"{self.name} ({self.symbol}): Mass={self.mass} MeV/c^2, Charge={self.charge}e, Spin={self.spin}, String Mode={self.string_mode}, Velocity={self.velocity}, Energy={self.energy} MeV, Momentum={self.momentum} MeV/c"

Velocity is a vector quantity representing the particle’s motion through space. To handle relativistic effects, the class includes methods to calculate the Lorentz factor, ensuring that velocities approach but do not exceed the speed of light (c). By integrating these properties, the Particle class forms the cornerstone of our simulation, enabling realistic modeling of particle behavior.

def lorentz_factor(self):
    vx, vy = self.velocity
    speed_squared = vx**2 + vy**2
    if self.mass == 0:
        return 1
    if speed_squared >= c**2:
        speed_squared = c**2 * 0.99
    return 1 / math.sqrt(1 - speed_squared / c**2)

def calculate_energy(self):
    vx, vy = self.velocity
    speed_squared = vx**2 + vy**2
    if self.mass == 0:
        return math.sqrt(speed_squared) * c
    gamma = self.lorentz_factor()
    return gamma * self.mass * c**2

def calculate_momentum(self):
    vx, vy = self.velocity
    speed_squared = vx**2 + vy**2
    if self.mass == 0:
        energy = math.sqrt(speed_squared) * c
        return (
            vx / math.sqrt(speed_squared) * energy / c,
            vy / math.sqrt(speed_squared) * energy / c,
        )
    gamma = self.lorentz_factor()
    return gamma * self.mass * vx, gamma * self.mass * vy

Modeling Fields

To simulate particle interactions accurately, we must consider the forces exerted by electric and magnetic fields. The ElectricField and MagneticField classes model these influences using the Lorentz force law.

The ElectricField class applies a constant electric field to particles, resulting in a force proportional to the particle’s charge and the field strength. This force induces acceleration, altering the particle’s velocity over time. The apply method updates the particle’s velocity vector, taking into account the time step of the simulation.

class ElectricField:
    def __init__(self, field_strength):
        self.field_strength = field_strength

    def apply(self, particle, time_step):
        force_x = particle.charge * self.field_strength[0]
        force_y = particle.charge * self.field_strength[1]

        if particle.mass > 0:
            ax = force_x / particle.mass
            ay = force_y / particle.mass

            new_vx = particle.velocity[0] + ax * time_step
            new_vy = particle.velocity[1] + ay * time_step
            particle.update_velocity(new_vx, new_vy)

Similarly, the MagneticField class models the effects of a magnetic field. The Lorentz force in this context is given by the cross product of the particle’s velocity and the magnetic field vector, scaled by the particle’s charge. This results in a perpendicular force that changes the particle’s trajectory without altering its speed. The apply method calculates this force and updates the particle’s velocity accordingly.

class MagneticField:
    def __init__(self, field_strength):
        self.field_strength = field_strength

    def apply(self, particle, time_step):
        vx, vy = particle.velocity
        Bz = self.field_strength[2]
        force_x = particle.charge * vy * Bz
        force_y = -particle.charge * vx * Bz

        if particle.mass > 0:
            ax = force_x / particle.mass
            ay = force_y / particle.mass

            new_vx = particle.velocity[0] + ax * time_step
            new_vy = particle.velocity[1] + ay * time_step
            particle.update_velocity(new_vx, new_vy)

Together, these classes provide a framework for simulating the dynamic interactions between particles and fields. By updating particle velocities and positions in response to these forces, we can capture the intricate dance of particles under the influence of electromagnetic fields, laying the foundation for more complex simulations.

Implementing Interactions – Decay, Collision, and Scattering

Particle Decay

Particle decay is a fundamental process in particle physics where unstable particles transform into other particles over time. This transformation is probabilistic in nature, governed by decay constants unique to each particle. In our simulation, the decay method models this behavior.

The method first checks if the particle has a non-zero decay constant. It then calculates the probability of decay occurring within the given time step using an exponential decay law. If decay occurs, the particle is replaced by its decay products, distributing the energy and momentum among them. This is crucial for accurately modeling particle lifetimes and the resultant particles from decays.

def decay(self, time_step):
    if self.decay_constant == 0:
        return [self]
    decay_probability = 1 - math.exp(-self.decay_constant * time_step)
    if random.random() < decay_probability:
        decay_products = random.choice(self.decays_to)
        total_energy = self.energy
        total_momentum = self.momentum
        num_products = len(decay_products)
        for product in decay_products:
            product.energy = total_energy / num_products
            product.momentum = (
                total_momentum[0] / num_products,
                total_momentum[1] / num_products,
            )
            if product.mass != 0:
                p_velocity_magnitude = math.sqrt(
                    (product.momentum[0] ** 2 + product.momentum[1] ** 2)
                    / product.mass**2
                )
                if p_velocity_magnitude >= c:
                    scale = (c * 0.99) / p_velocity_magnitude
                    p_velocity_magnitude *= scale
                product.update_velocity(
                    product.momentum[0] / product.mass,
                    product.momentum[1] / product.mass,
                )
            else:
                product.velocity = (
                    product.momentum[0] / product.energy * c,
                    product.momentum[1] / product.energy * c,
                )
        print(
            f"{self.name} decays into {', '.join([p.symbol for p in decay_products])}"
        )
        return decay_products
    else:
        return [self]

Collisions and Annihilation

Collisions between particles are another critical aspect of particle interactions. When two particles collide, they can form a new particle with combined mass and charge. The collide method handles this by creating a new particle and determining its properties based on the colliding particles.

Annihilation occurs when a particle meets its antiparticle, resulting in their mutual destruction and the creation of photons. This process is modeled in the annihilate method, which checks for charge neutrality and generates photon products accordingly.

def collide(self, other):
    new_mass = self.mass + other.mass
    new_charge = self.charge + other.charge
    new_velocity = (
        (self.velocity[0] + other.velocity[0]) / 2,
        (self.velocity[1] + other.velocity[1]) / 2,
    )
    new_particle = Particle(
        f"NewParticle({self.symbol}+{other.symbol})",
        f"{self.symbol}{other.symbol}",
        new_mass,
        new_charge,
        0,
        "mode_new",
        velocity=new_velocity,
    )
    print(
        f"{self.name} and {other.name} collide to form {new_particle.name}"
    )
    return new_particle

def annihilate(self, other):
    if self.charge + other.charge == 0:
        total_energy = self.energy + other.energy
        photon1 = Particle(
            "Photon", "γ", 0, 0, 1, "mode_γ", velocity=(c, 0)
        )
        photon2 = Particle(
            "Photon", "γ", 0, 0, 1, "mode_γ", velocity=(-c, 0)
        )
        photon1.energy = total_energy / 2
        photon2.energy = total_energy / 2
        print(
            f"{self.name} and {other.name} annihilate to form two photons"
        )
        return [photon1, photon2]
    else:
        return [self, other]

Scattering Events

Scattering events occur when particles interact and exchange velocities without forming new particles. This is common in high-energy physics where particles deflect off each other due to their electric charges. The scatter method models this by swapping the velocities of the interacting particles, thus simulating the scattering process.

def scatter(self, other):
    self_velocity = self.velocity
    other_velocity = other.velocity
    self.update_velocity(*other_velocity)
    other.update_velocity(*self_velocity)
    print(f"{self.name} and {other.name} scatter, exchanging velocities")

Integration of Interactions in the Simulation

To incorporate these interactions into our simulation, the simulate_interactions function checks for possible decays, collisions, and scattering at each time step. It updates the state of each particle accordingly, ensuring a realistic representation of particle dynamics.

def simulate_interactions(particles, fields, steps):
    for step in range(steps):
        print(f"\nStep {step + 1}")

        particles_copy = particles[:]
        interactions_done = set()
        for i in range(len(particles_copy)):
            for j in range(i + 1, len(particles_copy)):
                p1 = particles_copy[i]
                p2 = particles_copy[j]
                if (p1, p2) in interactions_done or (
                    p2,
                    p1,
                ) in interactions_done:
                    continue
                interactions_done.add((p1, p2))
                if p1.charge + p2.charge == 0 and random.random() < 0.1:
                    new_particles = p1.annihilate(p2)
                    particles.extend(new_particles)
                    if p1 in particles:
                        particles.remove(p1)
                    if p2 in particles:
                        particles.remove(p2)
                elif random.random() < 0.1:
                    new_particle = p1.collide(p2)
                    particles.append(new_particle)
                    if p1 in particles:
                        particles.remove(p1)
                    if p2 in particles:
                        particles.remove(p2)
                elif random.random() < 0.1:
                    p1.scatter(p2)

        for field in fields:
            for particle in particles:
                field.apply(particle, 1 / steps)

        for particle in particles[:]:
            new_vx = particle.velocity[0] + random.uniform(-0.01 * c, 0.01 * c)
            new_vy = particle.velocity[1] + random.uniform(-0.01 * c, 0.01 * c)
            particle.update_velocity(new_vx, new_vy)
            particle.update_position(1 / steps)
            decay_products = particle.decay(1 / steps)
            if decay_products != [particle]:
                particles.remove(particle)
                particles.extend(decay_products)

        print("Current particles:")
        for particle in particles:
            print(particle)

By accurately modeling decays, collisions, and scattering events, our simulation captures the complex interactions that define particle physics. This sets the stage for visualizing these dynamics and gaining deeper insights into the behavior of fundamental particles.

Simulating Particle Dynamics

Setting Up the Simulation

To effectively simulate particle interactions, we must initialize the simulation with a diverse set of particles and the fields that influence them. Our initial setup includes defining particles like electrons, protons, muons, and various neutrinos, each with specific properties such as mass, charge, and initial velocity. Additionally, we establish the electric and magnetic fields that will act upon these particles.

def simulate():
    electron = Particle(
        "Electron",
        "e",
        0.511,
        -1,
        1 / 2,
        "mode_1",
        velocity=(0.1 * c, 0.2 * c),
    )
    electron_neutrino = Particle(
        "Electron Neutrino",
        "νe",
        0.0000022,
        0,
        1 / 2,
        "mode_2",
        velocity=(0.1 * c, 0.2 * c),
    )
    muon_neutrino = Particle(
        "Muon Neutrino",
        "νμ",
        0.17,
        0,
        1 / 2,
        "mode_3",
        velocity=(0.1 * c, 0.2 * c),
    )
    tau_neutrino = Particle(
        "Tau Neutrino",
        "ντ",
        15.5,
        0,
        1 / 2,
        "mode_4",
        velocity=(0.1 * c, 0.2 * c),
    )
    muon = Particle(
        "Muon",
        "μ",
        105.66,
        -1,
        1 / 2,
        "mode_5",
        0.1,
        decays_to=[[electron, electron_neutrino, muon_neutrino]],
        velocity=(0.05 * c, 0.1 * c),
    )
    tau = Particle(
        "Tau",
        "τ",
        1776.86,
        -1,
        1 / 2,
        "mode_6",
        0.2,
        decays_to=[[muon, muon_neutrino, tau_neutrino]],
        velocity=(0.01 * c, 0.02 * c),
    )
    photon = Particle("Photon", "γ", 0, 0, 1, "mode_7", velocity=(c, c))
    proton = Particle(
        "Proton", "p", 938.27, 1, 1 / 2, "mode_8", velocity=(0.1 * c, 0.2 * c)
    )
    graviton = Particle(
        "Graviton", "G", 0, 0, 2, "mode_9", velocity=(0.5 * c, 0.5 * c)
    )
    axion = Particle(
        "Axion", "a", 1e-5, 0, 0, "mode_10", velocity=(0.05 * c, 0.1 * c)
    )
    particles = [electron, proton, tau, muon, photon, graviton, axion]
    fields = [ElectricField((1e5, 0)), MagneticField((0, 0, 1))]
    steps = 10

    print("Simulating interactions between particles in fields:")
    simulate_interactions(particles, fields, steps)
    visualize_particles(particles)

if __name__ == "__main__":  # pragma: no cover
    simulate()

Running the Simulation

The core of our simulation lies in the simulate_interactions function, which orchestrates the interactions between particles and their response to fields over multiple time steps. This function iterates through a predefined number of steps, updating particle states at each step.

During each iteration, the simulation checks for possible interactions between particles, such as decay, collision, and scattering. Fields are applied to each particle, updating their velocities according to the forces exerted. The particles' positions are also updated to reflect their movement through space.

def simulate_interactions(particles, fields, steps):
    for step in range(steps):
        print(f"\nStep {step + 1}")

        particles_copy = particles[:]
        interactions_done = set()
        for i in range(len(particles_copy)):
            for j in range(i + 1, len(particles_copy)):
                p1 = particles_copy[i]
                p2 = particles_copy[j]
                if (p1, p2) in interactions_done or (
                    p2,
                    p1,
                ) in interactions_done:
                    continue
                interactions_done.add((p1, p2))
                if p1.charge + p2.charge == 0 and random.random() < 0.1:
                    new_particles = p1.annihilate(p2)
                    particles.extend(new_particles)
                    if p1 in particles:
                        particles.remove(p1)
                    if p2 in particles:
                        particles.remove(p2)
                elif random.random() < 0.1:
                    new_particle = p1.collide(p2)
                    particles.append(new_particle)
                    if p1 in particles:
                        particles.remove(p1)
                    if p2 in particles:
                        particles.remove(p2)
                elif random.random() < 0.1:
                    p1.scatter(p2)

        for field in fields:
            for particle in particles:
                field.apply(particle, 1 / steps)

        for particle in particles[:]:
            new_vx = particle.velocity[0] + random.uniform(-0.01 * c, 0.01 * c)
            new_vy = particle.velocity[1] + random.uniform(-0.01 * c, 0.01 * c)
            particle.update_velocity(new_vx, new_vy)
            particle.update_position(1 / steps)
            decay_products = particle.decay(1 / steps)
            if decay_products != [particle]:
                particles.remove(particle)
                particles.extend(decay_products)

        print("Current particles:")
        for particle in particles:
            print(particle)

Monitoring Particle Evolution

Accurately capturing the evolution of each particle is critical to our simulation's fidelity. As particles move and interact, their velocities and positions are updated. These updates reflect changes due to forces from electric and magnetic fields, as well as interactions with other particles.

The update_velocity and update_position methods of the Particle class handle these updates. The update_velocity method ensures that particle speeds remain within relativistic limits, while the update_positionmethod tracks particle positions over time, enabling us to visualize their trajectories later.

def update_velocity(self, vx, vy):
    speed_squared = vx**2 + vy**2
    if speed_squared >= c**2:
        scale = math.sqrt(c**2 * 0.99 / speed_squared)
        vx *= scale
        vy *= scale
    self.velocity = (vx, vy)
    self.energy = self.calculate_energy()
    self.momentum = self.calculate_momentum()

def update_position(self, time_step):
    x, y = self.positions[-1]
    vx, vy = self.velocity
    new_x = x + vx * time_step
    new_y = y + vy * time_step
    self.positions.append((new_x, new_y))

By meticulously updating the state of each particle and accurately simulating their interactions, we lay the groundwork for visualizing particle dynamics. This approach allows us to observe the complex behaviors and interactions of particles as they move through space and interact with fields and each other.

Visualizing Particle Trajectories

Importance of Visualization in Particle Physics

Visualizing particle interactions is crucial for understanding the complex behaviors and dynamics that occur at the subatomic level. By observing the trajectories of particles, we can gain insights into how they move, interact, and change over time. Visualization transforms abstract data into a tangible form, making it easier to interpret and analyze the results of our simulations.

In particle physics, visualization helps bridge the gap between theoretical models and experimental observations. It allows researchers to validate their simulations against real-world data and provides a powerful tool for communicating findings within the scientific community and to the broader public.

Implementing Visualization

To visualize particle trajectories, we use the matplotlib library, a versatile tool for creating static, animated, and interactive plots in Python. Our approach involves plotting the positions of particles over time, resulting in a graphical representation of their paths through space.

We begin by defining a visualize_particles function, which takes a list of particles as input and generates a 2D plot of their trajectories. This function iterates through each particle's stored positions and plots them on a graph.

import matplotlib.pyplot as plt

def visualize_particles(particles):
    fig, ax = plt.subplots()
    for particle in particles:
        positions = zip(*particle.positions)
        ax.plot(*positions, label=particle.symbol)
    ax.set_xlabel("X Position (m)")
    ax.set_ylabel("Y Position (m)")
    ax.set_title("Particle Trajectories")
    ax.legend()
    plt.show()

In the visualize_particles function, we create a new figure and axis using plt.subplots(). We then iterate over the particles, extracting their positions and plotting them using ax.plot(). Each particle is labeled with its symbol for easy identification. Finally, we set the axis labels and title, and display the legend and plot.

Example Visualization

Let's walk through an example to illustrate how our visualization works. Consider a simulation with particles such as electrons, protons, and photons moving under the influence of electric and magnetic fields. After running the simulation, we call the visualize_particles function to generate a plot of their trajectories.

def simulate():
    electron = Particle(
        "Electron",
        "e",
        0.511,
        -1,
        1 / 2,
        "mode_1",
        velocity=(0.1 * c, 0.2 * c),
    )
    proton = Particle(
        "Proton", "p", 938.27, 1, 1 / 2, "mode_8", velocity=(0.1 * c, 0.2 * c)
    )
    photon = Particle("Photon", "γ", 0, 0, 1, "mode_7", velocity=(c, c))
    particles = [electron, proton, photon]
    fields = [ElectricField((1e5, 0)), MagneticField((0, 0, 1))]
    steps = 10

    print("Simulating interactions between particles in fields:")
    simulate_interactions(particles, fields, steps)
    visualize_particles(particles)

if __name__ == "__main__":  # pragma: no cover
    simulate()

Running this simulation, we can observe the paths of the electron, proton, and photon as they move through the fields and interact with each other. The resulting plot provides a clear visual representation of their trajectories, helping us understand how the electric and magnetic forces influence their motion.

By implementing and utilizing visualization techniques, we transform our simulation data into a format that is easy to analyze and interpret. This enables us to gain deeper insights into particle dynamics and enhances our ability to communicate our findings effectively.

Conclusion

Simulating and visualizing particle interactions offers invaluable insights into the fundamental processes that govern the universe. By leveraging Python to model particles and their interactions with electric and magnetic fields, we can observe the complex behaviors of particles, including decays, collisions, and scattering events. Visualization using matplotlibtransforms our simulations into tangible representations, making it easier to interpret and analyze the results. This approach not only enhances our understanding of particle physics but also demonstrates the power of computational tools in advancing scientific research. As we continue to refine our simulations, we move closer to uncovering the mysteries of the subatomic world.