Stochastic modelling of cosmic-ray sources for Galactic diffuse emissions

This paper employs stochastic Monte Carlo simulations to demonstrate that the discreteness of supernova remnant cosmic-ray sources introduces significant, energy-dependent uncertainties in Galactic diffuse gamma-ray and neutrino emission predictions, particularly in time-dependent diffusion scenarios, which may help reconcile theoretical models with high-energy observations from LHAASO and future experiments.

The Gist

Imagine our Galaxy, the Milky Way, as a giant, dark city. Invisible to the naked eye, this city is constantly being bombarded by high-speed particles called Cosmic Rays. These particles are like tiny, hyper-fast bullets zipping through the streets.

When these "bullets" crash into the gas and dust floating between the stars (the Interstellar Medium), they create a glow. This glow is Galactic Diffuse Emission (GDE). It's the "fog" of the galaxy, visible in gamma rays and neutrinos. Astronomers use this fog to figure out where the bullets are coming from and how fast they are moving.

For a long time, scientists modeled this fog as if the bullets were coming from a smooth, continuous sprinkler system spread evenly across the city. But in reality, the bullets don't come from a sprinkler; they come from specific, discrete "factories" (like exploding stars called Supernova Remnants).

This paper asks a simple but profound question: What happens to our map of the galaxy's fog if we realize the factories are individual, scattered points rather than a smooth blanket?

Here is the breakdown of their findings using some everyday analogies:

1. The "Smooth Fog" vs. The "Patchy Fog"

• The Old Way (Smooth Model): Imagine a fog machine that sprays mist evenly everywhere. If you look at the city, the fog looks the same everywhere, just thicker near the center. This is what previous models assumed.
• The New Way (Stochastic Model): Now, imagine the fog is actually made of steam rising from individual manholes scattered randomly across the city. If you stand near a manhole, it's thick and hot. If you stand between them, it's thin.
• The Finding: The authors ran a massive computer simulation (like running the city 1,000 times with the manholes in different random spots) to see how much the "fog map" changes.

2. Three Different "Factory Rules"

The scientists tested three different rules for how these factories (Supernova Remnants) release their particles:

• Scenario A: The "Burst" Factory

  • The Analogy: A factory explodes all its product at once, like a firework.
  • The Result: The fog has some patches, but they are relatively small. The difference between the "smooth map" and the "patchy map" is usually just a few percent. It's like noticing a few extra steam vents on a street, but the overall weather forecast doesn't change much.

• Scenario B: The "Smart" Factory

  • The Analogy: This factory releases heavy, slow particles first, and light, fast particles later.
  • The Result: Similar to the Burst scenario. The patches of fog are a bit different at different energies, but the overall picture remains fairly stable. The "smooth" model is still a pretty good guess.

• Scenario C: The "Trapped" Factory (Time-Dependent Diffusion)

  • The Analogy: Imagine the factory is inside a giant, sticky trap. The particles get stuck right next to the factory for a long time before they can escape into the city.
  • The Result: Chaos. The fog becomes incredibly patchy. Near a factory, the fog is thousands of percent thicker than the average. Far away, it's almost empty.
  • Why it matters: In this scenario, the "smooth map" is completely wrong. If we assume the fog is smooth, we miss the massive, localized storms of radiation right next to these factories.

3. The "LHAASO" Connection

We have new telescopes (like LHAASO) that can see this fog with incredible clarity, especially at very high energies (tens of Teravolts).

• The Problem: Sometimes, the smooth models predict less fog than what LHAASO actually sees. Scientists were worried they were missing something big in their physics.
• The Solution: The authors found that if the "Trapped Factory" scenario is true, the random "patches" of fog could explain why LHAASO sees more radiation than expected. The "extra" fog isn't a mystery; it's just the glow from a few nearby, super-active factories that the smooth model averaged out.

4. The "Correlation" Test

The scientists also checked if the fog looks the same at low energy (like a warm mist) and high energy (like a scalding steam).

• In the "Smooth" world: Low and high energy fog always look identical (just different colors).
• In the "Patchy" world:
  • In the "Burst" and "Smart" scenarios, the low and high energy fog still look mostly alike.
  • In the "Trapped" scenario, the low and high energy fog look completely different! The high-energy fog is clumped tightly around the factories, while the low-energy fog has spread out. This is a huge clue. If we see the fog looking different at different energies, it tells us the particles are getting "trapped" near their sources.

The Big Takeaway

This paper is a wake-up call for astronomers. We can no longer treat the galaxy's radiation fog as a smooth, uniform blanket.

• At low energies: The smooth model is a decent approximation.
• At high energies: The "patchiness" of individual sources becomes the most important factor.

The Analogy of the Future:
Imagine trying to predict the temperature of a room.

• Old View: You assume the heater is a giant radiator covering the whole wall. You predict a steady 70°F everywhere.
• New View: You realize the heater is actually a single, powerful blowtorch sitting in the corner.
• The Result: If you stand in the corner, you are burning up (100°F+). If you stand across the room, you are freezing (50°F). The "average" temperature of 70°F is useless for anyone actually living in that room.

Conclusion:
As our telescopes get sharper (higher spatial resolution), we will start seeing these "hot spots" of radiation. By studying them, we won't just be mapping the galaxy; we will be able to locate the specific "factories" (Supernova Remnants) that are accelerating these cosmic rays, and understand the physics of how they trap and release their energy. The randomness of the universe isn't just noise; it's a signal waiting to be decoded.

Technical Summary

Here is a detailed technical summary of the paper "Stochastic modelling of cosmic-ray sources for Galactic diffuse emissions" by Anton Stall and Philipp Mertsch.

1. Problem Statement

Galactic Diffuse Emissions (GDEs) in gamma rays and neutrinos arise from the interaction of Cosmic Rays (CRs) with the Interstellar Medium (ISM). While GDEs probe the CR intensity away from the Solar System, standard models typically assume a smooth, continuous distribution of CR sources (e.g., Supernova Remnants, SNRs).

However, CR sources are inherently discrete and localized. The exact positions and ages of these sources are unknown. This discreteness introduces a stochastic component to the CR intensity, particularly at high energies where escape timescales are short, and only a few nearby sources dominate the flux. The paper addresses the lack of quantitative understanding regarding how this source stochasticity affects GDE predictions, specifically:

• How do discrete sources alter the morphology of GDE sky maps compared to smooth models?
• What is the statistical distribution of GDE intensities (do they follow Gaussian or heavy-tailed distributions)?
• How does this stochasticity impact model uncertainties when comparing to high-precision data from experiments like LHAASO, IceCube, and future observatories (SWGO)?

2. Methodology

The authors developed a stochastic Monte Carlo framework to simulate hadronic GDEs.

A. Cosmic Ray Transport Model:

• Geometry: A simplified Galactic halo model with a disk ($z \approx 0$) and a halo height $H = 4$ kpc.
• Transport Equation: Solved for CR protons with rigidities $> 10$ GV, neglecting advection and inelastic collisions. The equation is $\partial_t \psi - \kappa \nabla^2 \psi = Q$.
• Green's Function: The solution for a single source is derived analytically using the method of mirror charges to satisfy free-escape boundary conditions at $z = \pm H$.
• Source Distribution: Sources are drawn from a spiral model with a Galactic bar. Source ages are uniformly distributed based on a supernova rate of $0.03 \text{ yr}^{-1}$.

B. Source Injection and Transport Scenarios:
The study compares three distinct scenarios to test the sensitivity of GDEs to source physics:

1. Burst-like: All rigidities escape simultaneously; time-independent diffusion.
2. Energy-Dependent Escape (EDE): High-rigidity particles escape earlier than low-rigidity ones (power-law dependence), but diffusion remains time-independent.
3. Time-Dependent Diffusion (TDD): All particles escape simultaneously, but diffusion is suppressed near the source ($\kappa_{start} \approx 10^{28} \text{ cm}^2\text{s}^{-1}$) for a fixed time ($t_{change} = 10$ kyr) before transitioning to standard rigidity-dependent diffusion. This mimics self-generated turbulence.

C. GDE Calculation:

• Hadronic Emission: Calculated via line-of-sight integration of CR proton intensity interacting with Galactic gas (HI and H2).
• Absorption: Includes photon-photon pair production ($\gamma\gamma \to e^+e^-$) using precomputed coefficients from GALPROP.
• Computational Efficiency: Utilized GPU acceleration (via jax) to handle the computationally expensive summation of millions of sources and matrix multiplications for line-of-sight integrals.

D. Statistical Approach:

• Monte Carlo Realizations: Simulated 1,000 full-sky realizations for morphology analysis and $10^6$ realizations for specific lines of sight (LOS) to study probability density functions (PDFs).
• Comparison: Contrasted discrete source results against a "smooth" model (the ensemble mean of the discrete realizations).

3. Key Contributions

1. Statistical Characterization of GDEs: Demonstrated that GDE intensity distributions are not purely Gaussian. They follow stable laws (Lévy distributions) with divergent variances in the Galactic disk, transitioning to mixtures of stable laws and Gaussians for LOS out of the disk.
2. Quantification of Model Uncertainty: Provided a rigorous quantification of the uncertainty introduced by source discreteness across different energy ranges and transport models.
3. Correlation Analysis: Investigated the correlation between low-energy (10 GeV) and high-energy (100 TeV) GDEs, showing how discrete sources break the strong correlation expected in smooth models.
4. LHAASO Reconciliation: Showed that source stochasticity in specific transport models (TDD) can bridge the gap between theoretical predictions and LHAASO measurements at high energies.

4. Key Results

A. Morphology and Sky Maps:

• Burst-like & EDE: Deviations from the smooth model are typically tens of percent. The morphology of excesses at 10 GeV and 100 TeV is correlated in the burst-like scenario but decoupled in the EDE scenario.
• Time-Dependent Diffusion (TDD): Shows extreme deviations, reaching order unity or larger (thousands of percent at 100 TeV). The confinement of CRs near sources creates distinct, localized "hotspots" that blur the line between diffuse and discrete emission.

B. Statistical Distributions:

• Stable Laws: In the Galactic plane, GDE intensities follow stable laws with stability indices $\alpha \approx 4/3$ (2D disk dominance) or $\alpha \approx 5/3$ (3D box dominance), depending on the intensity level.
• Mixture Models: For LOS out of the plane, the distribution is a mixture of a stable law (from nearby sources) and a Gaussian (from distant, numerous sources).
• Causality Effects: At very high energies, the implementation of causality (light-cone constraints) causes deviations from pure stable law predictions in the high-intensity tails.

C. Impact on Observables (Sky Windows & Profiles):

• Averaged Spectra (LHAASO Windows):
  • In Burst-like and EDE scenarios, the uncertainty due to source stochasticity is subdominant (<10%) compared to other model uncertainties (e.g., cross-sections, gas density). It cannot explain the discrepancy between smooth models and LHAASO data.
  • In the TDD scenario, stochasticity becomes a sizeable uncertainty (exceeding 50% at tens of TeV). The mean spectrum increases, and the uncertainty bands widen, potentially reconciling model predictions with LHAASO measurements in the inner Galaxy.
• Longitude Profiles: When spatial resolution is increased (e.g., longitude profiles), the influence of discrete sources becomes dominant. Individual realisations show significant structure (peaks) that smooth models miss.

D. Energy Correlations:

• Smooth Models: Predict a near-perfect correlation ($\rho \approx 1$) between 10 GeV and 100 TeV sky maps.
• Discrete Models:
  • Burst-like/EDE: The deviations ($\Delta J$) from the mean are poorly correlated between energies because different sources dominate at different energies.
  • TDD: The deviations are well-correlated because the initial transport phase is energy-independent, causing the same sources to imprint on both energy bands simultaneously. However, the total intensities ($J$) become poorly correlated due to the extreme enhancement of specific regions.

5. Significance and Conclusion

This work establishes that source stochasticity is a critical, non-negligible component of Galactic diffuse emission modeling, particularly at high energies and with high spatial resolution.

• For Data Interpretation: The "smooth" assumption is insufficient for interpreting high-energy data (LHAASO, SWGO) if the underlying transport involves time-dependent diffusion or if spatial resolution is high enough to resolve individual source imprints.
• Model Constraints: The lack of observed extreme "hotspots" in current data could constrain the rate of PeVatrons or rule out strong time-dependent diffusion scenarios. Conversely, the TDD scenario offers a mechanism to explain the excess flux observed by LHAASO without invoking new physics.
• Future Prospects: As experimental spatial resolution improves (especially beyond tens of TeV), measurements of GDEs will transition from constraining average source properties to locating individual CR sources and distinguishing between different injection/transport mechanisms.

The paper concludes that while stochasticity is subdominant for current low-resolution, broad-window analyses in standard models, it becomes the dominant source of uncertainty in time-dependent diffusion scenarios and for future high-resolution, high-energy observations.