Full-sky projection of diffuse gamma-ray sky.

Paper published: The diffuse gamma-ray sky of a Milky Way analog

Cosmic rays (CRs) are charged, relativistic particles, that constitute a significant energy component of the interstellar medium. They move relative to the thermal gas, which, combined with their inefficient cooling, allows them to establish long-lived pressure gradients that can have significant dynamical impacts on their host galaxies, for example by launching galactic outflows. They also affect the chemistry of gas through ionization and regulate star formation. The precise impact, however, is highly sensitive to the details of CR transport. This is a complex theoretical problem to untangle and currently largely unconstrained by observations.

Except in our local solar neighborhood, we observe CRs only indirectly through the radiation they produce, for example gamma-rays. The diffuse gamma-ray sky of our Milky Way is dominated by emission originating from neutral pions, which have decayed into gamma-ray photons. The pions themselves are created in collisions between CR protons and gas particles. Observations of the diffuse gamma-ray sky of the Milky Way provide an indirect probe of the underlying CR distribution, and traces the most dynamically important CRs.

To progress our understanding of CRs and constrain their impact from the simulation side, we need to produce synthetic observations and see how well they agree with real ones. Simulations also allow us to investigate the role of the local environment on the observed gamma-ray sky in a quantitative way, something that has not been attempted before in large-scale, self-consistent galaxy simulations.

For this project we ran a full simulation of an isolated Milky Way-like galaxy from the Rhea suite, including both magnetic fields and a relativistic CR fluid that is advected with and diffuses relative to the gas, where diffusion is anisotropic and directed along the magnetic field lines. The gamma-ray emission from CR protons is calculated in a post-processing step and assumes that the CRs are in a steady-state.

The left side shows a top-down view of the gas density in the galaxy plane. For six representative locations in the simulated Milky Way-like galaxy, the diffuse gamma-ray sky from CR protons is accurately modelled. The different realizations all show different unique features in the emission, pointing to the importance of the local environment.

We pick locations in our galaxy that mimic our own location in the Local Bubble, and compute the full diffuse gamma-ray sky as seen by a simulated observer at these points. The result is an incredibly diverse set of gamma-ray skies with unique filaments and features, especially at higher latitudes, see above (also Fig. 5 in the paper). This shows how important the local environment is for setting the features we see, which we quantify in detail in the paper.

The theoretically modelled gamma-ray sky spectrum is compared with observations of the gamma-ray sky as observed in the Milky Way. There is incredibly good aligment of the slope of the spectrum. See Fig. 10 in the paper. © The Authors (2026), CC-BY 4.0

One of our most surprising results is how well our simulation reproduces key observational properties of the Milky Way gamma-ray sky, both in terms of the total luminosity and flux and the shape of the gamma-ray spectrum, see above (and Fig. 10 in the paper). This is despite of the simplifying assumptions involved in our CR and gamma-ray modeling. For example, the CR transport parameters have been picked to match local CR observations, even though the underlying CR microphysics is almost certainly more complex. Our results suggest that our large-scale treatment captures a lot of the essential physics and provides an important benchmark for future studies with more sophisticated CR models. This is currently work in progress.

The diffuse gamma-ray sky of a Milky Way analog: Local diversity and global constraints. Karin Kjellgren, Philipp Girichidis, Maria Werhahn, Ralf S. Klessen, Christoph Pfrommer, Juan Soler, Brian Reville, Jim Hinton, Patrick Hennebelle, Noé Brucy and Simon C. O. Glover

A&A, 710 (2026) A163
DOI: https://doi.org/10.1051/0004-6361/202658922

© The Authors (2026). CC-BY 4.0

evo-step100-cutout

Properties of Isothermal Turbulence

In fluid dynamics, fluid motions can be separated into two very distinct types of flows. In Laminar flows, the fluid particles move in ordered Layers with little to no mixing between different Layers. A simple example of a laminar flow would be a viscous fluid flowing through a pipe. In contrast, turbulent flows are chaotic in nature. They can be observed in many aspects of everyday life, but are also of particular importance for many aspects of modern physics. Examples of this would be oceanic mixed layers, or chemical mixing in giant molecular clouds. While the statistical properties of turbulence are relatively well understood, mixing processes are still a field of active research, especially in astrophysical conditions.

There are two major simulation methods to model hydrodynamic systems. Many commonly used codes like AREPO use a mesh-based approach, where the simulation domain is subdivided into many smaller cells. Alternatively, smoothed particle hydrodynamics (SPH) based codes introduce particles that can move through the simulation domain and carry all necessary quantities like pressure and density without requiring any grid. Such SPH-based codes are particularly well suited to study mixing problems. As the simulation method in Lagrangian in nature, mixing can simply be studied by tracing the SPH-particles without requiring any additional tracer particles.

For this reason, our group is interested in using the new SPH-EXA simulation code to study mixing problems related to giant molecular clouds and star formation. SPH-EXA is a new, GPU-accelerated simulation code that is developed by research teams at the universities in Basel and Zürich. It is specifically designed for high-resolution simulations and can handle billions of SPH-particles with relative ease. To achieve the necessary computational power, we use the Alps supercomputer, one of the largest supercomputers in the world.

So far, our work has mainly focused on validating and improving this new simulation code. This is especially important for subsonic turbulence, a regime with which the SPH method has historically struggled. Through several improvements in the simulation method, SPH-EXA can now produce accurate probability density functions, power spectra and structure functions in the regime of subsonic turbulence using the SPH method. Our results match those from established mesh-based codes like AREPO (see Cabézon et al (2025) for more details). This is a very important milestone, as we are now confident that our code accurately models turbulent flows in both the sub- and supersonic regime.

Slices of a turbulence simulation. The six panels show the density, velocity, divergence of velocity, curl of velocity and backward and forward finite-time Lyapunov Exponent.
Slices of a turbulence simulation. The six panels show the density, velocity, divergence of velocity, curl of velocity and backward and forward finite-time Lyapunov Exponent.

The SPH method is not only beneficial for studying mixing coefficients. It also gives easy access to diagnostics like the forward and backward in time finite-time Lyapunov exponent (FTLE). The FTLE is, broadly speaking, a measurment of chaos. Here, we use this quantity to identify the most attracting and repelling surfaces in a flow and is useful to visualize coherent structures in complex flows. The Figure above shows a slice through a simulation of isothermal, supersonic turbulence with periodic boundary conditions. This simulation was carried out with a Mach number of 4 and a resolution of 8 10 9 SPH-particles. One can see that both the forward and backward in time the high-FTLE regions correspond to the shock fronts in the fluid. This is expected, as the shock fronts are regions where the gas is compressed and then decompressed very rapidly.

With all this in mind, we can now move forward to studying mixing processes under astrophysical conditions. An specific application is mixing of hot and cold gas phases in the interstellar medium where cold gas is able to form stars but stars then heat gas due to their radiation and, finally, during their supernova explosion. These idealised turbulence simulations help to better understand how the different gas phases mix and interact to regulate star formation.

This work has been done during the Master Thesis of Oliver Avril.

A photograph of the LUMI-G supercomputer

European Supercomputer Aids Heidelberg Astrophysics

A Swiss-German research team hopes to unlock the secrets of star formation using Europe’s fastest computer – the LUMI-G supercomputer in Kajaani (Finland) run by an international consortium. The researchers have received computing time on the LUMI-G for the required simulations, which will also utilise a newly developed simulation code. Also Heidelberg astrophysicists are contributors to the research. Project partner Prof. Dr Ralf Klessen of Heidelberg University’s Centre for Astronomy (ZAH) also anticipates groundbreaking insights for his own research.