🔭 astrophysics

Implementation of CR Energy SPectrum (CRESP) algorithm in PIERNIK MHD code. II. Propagation of Primary and Secondary nuclei in a magneto-hydrodynamical environment

This paper presents an extension of the PIERNIK MHD code's CRESP module to simulate the coupled production and propagation of primary and secondary cosmic ray nuclei, revealing that increased diffusion coefficients and supernova rates respectively enhance and suppress the secondary-to-primary flux ratio due to magnetic buoyancy effects and reduced particle residence times.

Original authors: Antoine Baldacchino-Jordan, Michał Hanasz, Mateusz Ogrodnik, Dominik Wóltański, Artur Gawryszczak, Andrew W. Strong, Philipp Girichidis

Published 2026-03-16

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Antoine Baldacchino-Jordan, Michał Hanasz, Mateusz Ogrodnik, Dominik Wóltański, Artur Gawryszczak, Andrew W. Strong, Philipp Girichidis

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the universe not as a silent void, but as a bustling, chaotic city where invisible, high-speed particles called Cosmic Rays zip around like hyperactive commuters. These aren't just random travelers; they are the "ghosts" of exploded stars, carrying energy that shapes the very air (the Interstellar Medium) they fly through.

This paper is about building a super-sophisticated traffic simulator for these cosmic commuters. The authors, using a powerful computer code called PIERNIK, have upgraded their model to track not just the main commuters (primary particles like Carbon and Oxygen), but also the "accidental" passengers they create along the way (secondary particles like Lithium, Beryllium, and Boron).

Here is the breakdown of their discovery, explained through everyday analogies:

1. The Setup: A Cosmic City in a Box

The researchers built a digital "box" representing a slice of our galaxy near the Sun. Inside this box, they simulated:

The Gas: The air of the city.
The Magnetic Fields: Invisible train tracks guiding the particles.
The Supernovae: Random explosions (like fireworks) that inject new cosmic rays into the system.

In previous versions of their code, they could only track the "main" particles. In this new version, they added the ability to track what happens when a fast-moving primary particle smashes into a gas atom. It's like watching a billiard ball hit a cluster of balls; the original ball slows down, and new, smaller balls (secondary particles) scatter in different directions.

2. The Big Surprise: The "Leaky Box" vs. The "Living City"

For decades, scientists used a simple model called the "Leaky Box." Imagine a bucket with a hole in the bottom.

Old Logic: If you make the hole bigger (faster diffusion), particles escape faster. Therefore, you have fewer particles inside, and fewer collisions, meaning fewer "accidental" secondary particles.
The New Discovery: The authors found that in their living, breathing simulation, the opposite happens!
- The Analogy: Imagine the cosmic rays are hot air balloons. When they move too fast (high diffusion), they don't just escape; they actually push the air around them. This creates a "buoyancy" effect that changes the shape of the magnetic "train tracks."
- The Result: When the diffusion is high, the magnetic tracks become more horizontal and less vertical. This traps the primary particles near the "ground" (the galactic disk) for longer because they can't easily climb out.
- The Consequence: Because the primary particles stay in the crowded disk longer, they smash into gas atoms more often, creating more secondary particles. So, faster diffusion = more secondary particles. This flips the old theory on its head!

3. The Supernova Factor: The "Wind" Effect

The team also tested what happens if you change the number of explosions (Supernovae).

More Explosions: This creates a stronger "wind" blowing upward out of the galaxy. It's like a giant fan turning on.
The Effect: This wind sweeps the cosmic rays out of the disk very quickly. They don't stay long enough to crash into gas atoms and create secondaries.
The Result: More explosions actually lead to fewer secondary particles (a lower Boron-to-Carbon ratio) because the "commuters" are kicked out of the city too fast to cause accidents.

4. Why This Matters

The scientists compared their simulation results with real data collected by the AMS-02 instrument on the International Space Station.

They found that their complex, self-consistent model (where the particles and the gas talk to each other) behaves differently than the old, static models.
They discovered that to match the real universe, you can't just tweak one number (like the diffusion rate). You have to understand how the wind, the magnetic fields, and the explosions all dance together.

The Takeaway

This paper is a major step forward because it stops treating the galaxy as a static room and starts treating it as a dynamic ecosystem.

Old View: Cosmic rays are just dust blowing in the wind.
New View: Cosmic rays are the wind itself. They push the gas, bend the magnetic fields, and change their own path by doing so.

By understanding this complex dance, scientists can finally figure out exactly how long these particles stay in our galaxy and where they come from, using the "accidental" particles (like Boron) as clues to solve the mystery.

1. Problem Statement

Cosmic Ray (CR) propagation is traditionally modeled using phenomenological "leaky-box" models or stationary diffusion-advection equations (e.g., GALPROP) that assume a static Interstellar Medium (ISM). However, these models often neglect the dynamic feedback of CRs on the ISM, such as CR-driven winds and magnetic field restructuring. While self-consistent Magnetohydrodynamic (MHD) models exist, they often treat CRs as a spectrally unresolved "grey" fluid or focus only on electrons.

There is a need for a self-consistent, spectrally resolved model that simultaneously tracks the propagation of multiple primary and secondary CR nuclei species while dynamically coupling them to the evolving thermal gas, magnetic fields, and supernova (SN) feedback. Specifically, the authors aim to validate the production and propagation of secondary nuclei (Li, Be, B) resulting from the spallation of primary nuclei (C, N, O) within a realistic, dynamic galactic environment.

2. Methodology

Numerical Framework: PIERNIK MHD Code

The study utilizes PIERNIK, a grid-based, multifluid MHD code. The authors extend the existing CRESP (Cosmic Ray Energy SPectrum) algorithm, previously used for electrons, to include hadronic CR species.

Two-Moment Approach: The code solves the Fokker-Planck equation for CR distribution functions using a two-moment approach (number density $n_a$ and kinetic energy density $e_a$ ) within momentum bins.
Piece-wise Power-Law Approximation: Within each momentum bin, the distribution function is approximated as a power law. This allows for efficient numerical integration of energy gain/loss terms without resolving the full distribution function at every step.
Spectral Resolution: The model tracks 16 momentum bins ranging from non-relativistic to ultra-relativistic energies ( $5 \times 10^{-1}$ to $10^5$ proton masses).

Physical Processes Implemented

Primary Nuclei: $^{12}\text{C}$ , $^{14}\text{N}$ , and $^{16}\text{O}$ are injected via supernova remnants (SNRs). Protons are treated as a "grey" (spectrally unresolved) relativistic fluid to drive the ISM dynamics.
Secondary Nuclei: Secondary species ( $^{7}\text{Li}$ , $^{9}\text{Be}$ , $^{10}\text{Be}$ , $^{10}\text{B}$ , $^{11}\text{B}$ ) are generated via spallation (hadronic collisions) of primary nuclei against thermal ISM protons. Cross-sections are taken from GALPROP.
Radioactive Decay: The decay of unstable isotopes (specifically $^{10}\text{Be}$ ) is included.
Energy Losses: The model includes adiabatic expansion losses and Coulomb cooling (using the Bethe-Block formula). Streaming effects are currently omitted.
Transport Mechanisms: CRs undergo advection with the thermal gas and anisotropic, rigidity-dependent diffusion ( $D_{\parallel} \propto R^{\delta}$ ) along magnetic field lines.

Simulation Setup

Domain: A 3D gravity-stratified box ( $0.5 \times 1 \times 8$ kpc) representing the local ISM near the Solar System.
Initial Conditions: Hydrostatic equilibrium with $n_{\text{ISM}} \approx 1 \text{ cm}^{-3}$ and $T \approx 7000$ K.
Feedback: Randomly distributed supernovae inject CR protons (10% of $10^{51}$ erg per SN) to drive the system.
Parameter Study: The authors performed simulations varying three key parameters:
1. SN Rate: 20, 60, and 80 kpc $^{-2}$ Myr $^{-1}$ .
2. Diffusion Coefficient ( $D_{\parallel}$ ): $3 \times 10^{28}$ and $6 \times 10^{28}$ cm $^2$ s $^{-1}$ .
3. Rigidity Index ( $\delta$ ): 0.3 and 0.5.

3. Key Contributions

Extension of CRESP to Hadronic Species: This is the first implementation of the CRESP algorithm in PIERNIK to handle multiple spectrally resolved primary and secondary nuclei simultaneously, including spallation source terms and radioactive decay.
Self-Consistent CR-ISM Coupling: Unlike stationary models, this study demonstrates how CRs dynamically shape the ISM (via buoyancy and winds) and how that reshaped environment, in turn, alters CR transport.
Novel Discovery on Diffusion and B/C Ratio: The paper identifies a counter-intuitive relationship between the diffusion coefficient and the Boron-to-Carbon (B/C) ratio, driven by magnetic field restructuring.

4. Key Results

Dynamical Evolution

CR-Driven Outflows: CR injection drives vertical gas outflows and creates Parker instability structures (magnetic loops). Higher SN rates lead to stronger winds and faster vertical gas velocities.
Magnetic Field Restructuring: CR buoyancy creates vertical magnetic field components. The strength of these vertical fields depends on the diffusion efficiency.

Spectral Properties

Slope Deviations: The evolved CR spectra show deviations from the initial injection slope ( $q \approx 4.1$ $q \approx 4.1$ ).
- At relativistic energies, the spectral slope difference between primaries and secondaries ( $\Delta a_R$ ) matches the diffusion power index $\delta$ .
- At non-relativistic energies, slopes are steeper due to Coulomb cooling and spallation, with the disc showing harder spectra than the halo.
Energy Density Peaks: Spectra exhibit peaks around 10–30 GeV/nucleon due to the interplay of injection, cooling, and diffusion.

Observables: B/C and $^{10}\text{Be}/^{9}\text{Be}$ Ratios

Impact of SN Rate:
- Higher SN Rate $\rightarrow$ Lower B/C. Stronger winds increase advection, reducing the residence time of primary CRs in the dense disk where spallation occurs.
- Higher SN rates also lead to lower $^{10}\text{Be}/^{9}\text{Be}$ ratios due to faster escape and decay.
Impact of Diffusion Coefficient ( $D_{\parallel}$ ):
- Novel Finding: Contrary to standard leaky-box predictions (where higher $D$ implies lower B/C), this model shows higher $D_{\parallel}$ leads to higher B/C.
- Mechanism: Higher diffusion distributes CRs more evenly along horizontal field lines, reducing the local CR pressure gradients that drive magnetic buoyancy. This results in weaker vertical magnetic fields ( $\langle |B_z|/|\vec{B}| \rangle$ decreases). Weaker vertical fields trap primary CRs longer in the disk, increasing the production of secondaries (Boron).
Impact of Rigidity Index ( $\delta$ ):
- Increasing $\delta$ from 0.3 to 0.5 significantly steepens the B/C ratio at relativistic energies, bringing the simulation results into better agreement with AMS-02 observational data (which favors $\delta \approx 1/3$ ).

5. Significance and Limitations

Significance:

The study validates that self-consistent MHD simulations are necessary to accurately interpret CR observables. The "effective" diffusion experienced by CRs is not just a function of the input parameter $D_{\parallel}$ but is modulated by the magnetic field geometry, which is itself a product of CR feedback.
It provides a robust numerical framework for comparing complex CR propagation models with high-precision data from AMS-02 and Voyager.

Limitations & Future Work:

Resolution: The current grid resolution (62.5 pc) is coarser than state-of-the-art simulations (1–4 pc), potentially missing small-scale turbulence.
Simplified Physics: The model omits CR streaming, multi-phase ISM thermal structure, and detailed SN thermal/kinetic feedback.
Geometry: The simulation uses a local patch with periodic boundaries, which may introduce biases compared to a full galactic halo model.

Conclusion:
The authors successfully demonstrated that the PIERNIK MHD code, extended with the CRESP algorithm, can self-consistently model the production and propagation of CR nuclei. The results highlight that CR transport observables (like B/C) are sensitive to the dynamic coupling between CRs, gas, and magnetic fields, challenging the assumptions of static, phenomenological models.