Probabilistic modeling of Cherenkov emission from particle showers

This paper presents a probabilistic modeling approach using parameter distributions to efficiently simulate Cherenkov light yield fluctuations from particle showers in ice and water, offering a more accurate alternative to full Monte Carlo simulations for neutrino telescope experiments.

The Gist

The Big Picture: Predicting the "Light Show" of Cosmic Particles

Imagine you are trying to understand a massive fireworks display happening deep underwater. You can't see the fireworks directly, but you can see the ripples of light they create in the water. In the world of physics, these "fireworks" are particle showers.

When a high-energy particle (like a neutrino from deep space) smashes into an atom in ice or water, it doesn't just stop. It shatters the atom and creates a cascade of hundreds of smaller, lower-energy particles. These particles zoom through the ice faster than light can travel in that ice, creating a cone of blue light called Cherenkov radiation. This is how telescopes like IceCube "see" the universe.

The Problem:
To understand what happened, scientists need to know exactly how that light is distributed. Is it a bright flash right at the start? Does it fade out slowly? Does it have a second burst later on?

Traditionally, scientists used a "blurry average." They would say, "On average, a 1 TeV particle makes a light show that looks like this." But in reality, every single particle shower is unique. Some are messy, some are clean, some have extra bursts of light, and some are dimmer than expected. Using the "average" is like trying to predict the weather by only looking at the monthly average temperature; it misses the storms and the heatwaves.

The Solution:
This paper introduces a new, smarter way to simulate these light shows. Instead of guessing the average, they built a probabilistic model. Think of it as moving from a "one-size-fits-all" t-shirt to a tailor-made suit for every single event.

---

How They Did It: The "Digital Lab"

1. The Heavy Lifting (FLUKA)
First, the team used a super-complex computer program called FLUKA. Imagine this as a hyper-realistic video game engine for physics. They ran millions of simulations, smashing particles into digital ice blocks to see exactly how the light was created in every single scenario.

• The Catch: Running these simulations is incredibly slow and expensive, like trying to render a 4K movie frame-by-frame for every single event. You can't do this for every real-world observation.

2. The Shortcut (The New Model)
Since they couldn't run the slow simulations for every future event, they took the data from their "heavy lifting" and built a statistical recipe.

• They realized that while every shower is unique, they all follow certain patterns.
• They created a model that doesn't just give you one answer, but a range of possible answers with specific odds.
• The Analogy: Instead of saying "The light will be 50 units bright," the model says, "There's a 70% chance it's 45-55 units, a 20% chance it's 30 units, and a 10% chance it's 60 units."

The Two Main Ingredients

The model breaks the light show down into two parts:

1. The "Volume" of Light (Amplitude)
How much total light is there?

• Old Way: Assumed the amount of light was a perfect bell curve (Gaussian).
• New Way: Discovered that the light amount is often "skewed." Sometimes, a particle decays early and creates almost no light; other times, it creates a massive burst. The new model uses a special mathematical shape (called a Skew Normal or Normal-Inverse Gaussian) to capture these weird, lopsided tails.
• Metaphor: Imagine a bag of marbles. The old model assumed the bag always had exactly 100 marbles. The new model knows that sometimes the bag has 90, sometimes 120, and sometimes it's empty because the marbles fell out early.

2. The "Shape" of the Light (Profile)
How is that light spread out along the path?

• Old Way: Assumed the light always fades out in a smooth, predictable curve.
• New Way: Used a technique called B-splines. Imagine a flexible ruler that can bend to fit any shape. The model learns how this "ruler" bends for different types of particles and energies.
• Metaphor: The old model drew a smooth, straight line for a road. The new model draws a road that can have potholes, bumps, and sharp turns, because that's what real particle showers actually look like.

Why Does This Matter?

1. Better Detective Work
Neutrino telescopes are trying to figure out where cosmic particles came from and what they are. If you use the "blurry average" model, you might misjudge the direction or energy of the particle. By using this new, fluctuating model, scientists can reconstruct the event much more accurately.

• Analogy: If you are trying to find a lost hiker by looking at their footprints, the old model assumes everyone leaves perfect, identical footprints. The new model knows that muddy boots, running, or slipping leave different, messy prints, helping you track the hiker more accurately.

2. Speed vs. Accuracy
The new model is a "best of both worlds" compromise. It is almost as accurate as the slow, heavy simulations but runs thousands of times faster. This allows scientists to simulate billions of events quickly, which is crucial for training Artificial Intelligence (AI) to recognize real signals from background noise.

The Bottom Line

This paper is about upgrading the toolkit for cosmic detectives. They realized that nature is messy and unpredictable. Instead of forcing nature into a neat, average box, they built a flexible, probabilistic model that embraces the chaos. This means that when the next big cosmic event happens, our telescopes will be better equipped to understand exactly what it was, where it came from, and what it tells us about the universe.

Technical Summary

1. Problem Statement

High-energy particle showers, initiated by cosmic rays or neutrino interactions in media like ice or water, are central to astroparticle physics. While the development of these showers is governed by the Standard Model and exhibits universality, accurately simulating them for neutrino telescopes (e.g., IceCube) presents significant challenges:

• Computational Cost: Full Monte Carlo (MC) simulations (e.g., GEANT4, FLUKA) that track every secondary particle are computationally prohibitive when scaled to the massive event rates required for next-generation experiments.
• Limitations of Current Approximations: Existing fast simulation methods rely on parametrized average profiles (typically using a single Gamma distribution). These deterministic models fail to capture event-to-event fluctuations in both the total Cherenkov light yield ($\hat{\ell}_{tot}$) and the longitudinal shower profile ($d\hat{\ell}/dx$).
• Impact on Physics: Ignoring these fluctuations leads to inaccuracies in reconstructing neutrino properties (flavor, energy, direction) and distinguishing signal from background, particularly for rare events or those with small separation distances.

2. Methodology

The authors propose a probabilistic framework that replaces deterministic averages with distributions of parameters, allowing for the generation of realistic, fluctuating shower profiles.

A. High-Fidelity Simulation (FLUKA)

• Tool: The study utilizes FLUKA (version 2025.1.0), a high-precision MC transport code, to generate a large dataset of particle showers in ice ($\rho_0 = 0.9216$ g/cm$^3$).
• Inputs: Simulations cover primary particles relevant to neutrino Deep Inelastic Scattering (DIS): electrons ($e^-$), photons ($\gamma$), and a wide range of hadrons (protons, neutrons, pions, kaons, hyperons, etc.) across energies from 1 GeV to 1 PeV.
• Physics Settings: The configuration enables advanced physics processes often omitted in standard parameterizations, including photonuclear and electronuclear interactions, heavy fragment evaporation, and coalescence.
• Output: The simulation calculates the weighted track length ($\hat{\ell}$) for charged particles, which serves as a proxy for Cherenkov light yield. The output is binned longitudinally to produce $d\hat{\ell}/dx$ profiles.

B. Probabilistic Modeling Strategy

The authors decompose the shower description into two independent components: Amplitude and Shape.

1. Shape Modeling ($a, b$ parameters):

   • Individual shower profiles are fitted to a Gamma distribution: $d\hat{\ell}/dx \propto x^{a-1}e^{-bx}$.
   • To handle the unbounded nature of parameters $a$ and $b$, they are transformed into a unit square using $a' = 1/\sqrt{a}$ and $b' = 1/(1+b^2)$.
   • B-Spline Density Estimation: A joint probability density function (PDF), $f(a', b'; E)$, is constructed using tensor products of penalized B-splines (a Generalized Linear Model approach). This captures the complex, energy-dependent correlations between shape parameters.
   • Smoothing: A penalty term is applied to prevent overfitting, optimized via the Akaike Information Criterion (AIC).

2. Amplitude Modeling ($\hat{\ell}_{tot}$):

   • The total weighted track length $\hat{\ell}_{tot}$ is modeled separately.
   • Electromagnetic (EM) Showers: Modeled using a Normal-Inverse Gaussian (NIG) distribution to account for observed skewness caused by nuclear effects.
   • Hadronic Showers: Modeled using a Skew Normal distribution.
   • Energy Dependence: The parameters of these distributions are fitted as functions of primary energy using 6th-degree polynomials in $\log_{10}(E)$.

3. Sampling:

   • New shower profiles are generated by sampling $\hat{\ell}_{tot}$ from the parametric distribution and $(a', b')$ from the B-spline PDF, then reconstructing the Gamma profile.

3. Key Contributions

• Event-to-Event Fluctuation Modeling: The paper moves beyond "average" shower profiles, providing a statistical method to generate unique, fluctuating shower shapes and amplitudes for every simulated event.
• Advanced Nuclear Physics Integration: By utilizing FLUKA with specific nuclear interaction settings (photonuclear/electronuclear), the model captures physical effects (e.g., nucleon emission) that cause significant deviations in EM shower yields, which were previously unmodeled in standard parameterizations.
• Hybrid Parametric Approach: The use of B-splines for shape parameters allows for a flexible, non-parametric representation of complex correlations without the computational cost of full MC tracking.
• Open-Source Implementation: The authors provide a Python package containing the fitted B-spline coefficients and polynomial parameters, making the model immediately usable for other experiments.

4. Results

• Accuracy Validation:
  • Wasserstein Distance ($W_1$): The new model significantly reduces the distance between simulated and modeled distributions compared to previous average-based parametrizations (Refs. [2, 3]).
  • Kolmogorov-Smirnov (KS) Test: Samples drawn from the model are statistically consistent with the FLUKA simulation data for shape parameters ($a', b'$) and amplitude ($\hat{\ell}_{tot}$) across most energies.
• Fluctuation Characteristics:
  • EM Showers: Exhibit lower fluctuations but show a distinct left-skew in $\hat{\ell}_{tot}$ due to nuclear effects, which the NIG model captures accurately.
  • Hadronic Showers: Show larger variance in both shape and amplitude. The model successfully reproduces the broadening of distributions and the "tail" effects seen in FLUKA.
  • Outliers: The model handles the majority of showers well. However, at low energies (<100 GeV for hadrons), rare decay channels (e.g., $K^0_S \to \pi^0\pi^0$) create multi-peak structures or zero-yield tails that the single-Gamma approximation cannot fully capture. The authors recommend a minimum energy threshold (100 GeV for hadrons, 10 GeV for EM) for this parametrization.
• Application to Neutrino Interactions:
  • Applied to 100 TeV $\nu_e$ and $\nu_\tau$ DIS interactions using PYTHIA8 for the initial interaction.
  • The model successfully generates complex, multi-peaked longitudinal profiles resulting from the superposition of multiple sub-showers (e.g., from $\tau$ decays or NC interactions), demonstrating its utility for reconstructing complex topologies.

5. Significance

• Improved Reconstruction: By incorporating realistic fluctuations, neutrino telescopes can achieve more accurate energy and direction reconstruction. Fluctuations in shower elongation can be exploited to better distinguish between shower types (e.g., EM vs. Hadronic).
• Background Rejection: Accurate modeling of rare, fluctuating shower shapes is crucial for identifying rare signals (like Double Bang events or specific BSM signatures) that might be obscured by incorrect background modeling.
• Computational Efficiency: The probabilistic model offers a "fast simulation" alternative that is orders of magnitude faster than full MC tracking while retaining the statistical fidelity of the underlying physics, enabling large-scale simulations for future experiments (e.g., IceCube-Gen2).
• Future-Proofing: The framework is designed to be easily updated with new physics settings or primary particle types, serving as a robust foundation for the next generation of neutrino astronomy.