Cosmic ray transport and acceleration in an evolving shock landscape

TL;DR

This paper advances our understanding of cosmic-ray transport and re-acceleration by modeling time-dependent shocks in the Galactic halo, motivated by local outbursts near the Galaxy's edge. It extends the CRPropa framework to include time-dependent winds and couples it to Athena++ hydrodynamics to capture both planar and Sedov–Taylor-type blast waves, deriving SDE-based transport and momentum evolution for CRs. The results show that re-accelerated CRs can propagate back to the Galaxy and contribute to the spectrum near and above the knee, with spectra that deviate from simple power laws; in particular, single shocks tend to yield flat-to-steepening spectra, while colliding shocks can produce transient hardening near collision times. Overall, the work highlights the importance of local, time-dependent shock landscapes for shaping the Galactic CR spectrum and demonstrates a versatile, validated methodology for exploring such scenarios with realistic hydrodynamics and diffusion physics.

Abstract

The sources of cosmic rays between the knee and the ankle are still debated. The Galactic wind and its termination shock have been proposed to contribute to this transition between Galactic and extragalactic origin, but another possibility is large-scale shock structures from local sources in the Milky Way. In this paper, we investigate CR transport in a time-dependent landscape of shocks in the Galactic halo. These shocks could result from local outbursts, e.g. starforming regions and superbubbles. CRs re-accelerated at such shocks can reach energies above the knee. Since the shocks are closer to the Galaxy than a termination shock and CRs escape downstream, they can propagate back more easily. With such outbursts happening frequently, shocks will interact. This interaction could adjust the CR spectrum, particularly for the particles that are able to be accelerated at two shocks simultaneously. The transport and acceleration of CRs at the shock is modeled by Stochastic Differential Equations (SDEs) within the public CR propagation framework CRPropa. We developed extensions for time-dependent wind profiles and for the first time connected the code to hydrodynamic simulations, which were run with the public Athena++ code. We find that, depending on the concrete realization of the diffusion tensor, a significant fraction of CRs can make it back to the Galaxy. These could contribute to the observed spectrum around and above the CR knee ($E \gtrsim 10\,\mathrm{PeV}$). In contrast to simplified models, a simple power-law does not describe the energy spectra well. Instead, for single shocks, we find a flat spectrum ($E^{-2}$) at low energies, which steepens gradually until it reaches an exponential decline. When shocks collide, the energy spectra transiently become harder than $E^{-2}$ at high energies.

Figures (21)

• Figure 1: Shock profiles in lab frame (left) and frame in which shock is stationary (right). Assuming a compression ratio $q = 4$ and that the preshock medium does not move ($v_1 = 0$), the postshock speed is $3/4$ of the shock speed. In the shock frame, upstream and downstream side can be specified and the downstream speed equals $1/q$ the upstream speed.
• Figure 2: Left: Wind profiles and space-energy histogram in lab frame. Pseudo-particles are injected in the undisturbed medium, the shock is moving through it. Right: Stationary wind profile and space-energy histogram. Pseudo-particles are injected at the shock. For both columns, the top row shows spectra and wind profiles at $T=50$ and the bottom row shows them at $T=90$.
• Figure 3: Energy spectra at the shocks in lab frame (moving with the shock) (blue) and in stationary shock frame (orange). The resulting spectral slope and time evolution are the same. Error bars are obtained from the Monte-Carlo error in each bin and proportional to $\sqrt{N}$ where $N$ is the number of independent candidates in the respective bin. The discrepancy mainly is due to the time integration in case of the stationary profile.
• Figure 4: Top: Wind profile of the self-similar Sedov-Taylor solution. The shock slows down over time $\propto t^{-3/5}$. Here $E_{\mathrm{inj}} = 20$, $\rho_0 = 1$. Athena++ simulation (solid line) is compared to analytical profile with different shock widths, $l_{\mathrm{sh}} = [0.1, 0.05, 0.01]$ (dashed, dash-dotted, dotted line). Bottom: Relative error between the Athena++ simulation and approximated analytical profile. For the smallest shock width the error increases over time since the shock position of the analytical solution and Athena++ simulation diverge.
• Figure 5: Acceleration at Sedov-Taylor blast wave. Particles are accelerated at the shock and are cooled in the expanding wind.