# Cosmic Rays in Superbubbles

T. Tolksdorf, I.A. Grenier, T. Joubaud, R. Schlickeiser

► To cite this version:

T. Tolksdorf, I.A. Grenier, T. Joubaud, R. Schlickeiser. Cosmic Rays in Superbubbles. The Astrophysical Journal, 2019, 879 (2), pp.66. 10.3847/1538-4357/ab24c6. http://arXiv.org/abs/1902.04662

HAL Id: hal-02223028

https://hal.science/hal-02223028v1

Submitted on 10 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Cosmic Rays in Superbubbles

T. Tolksdorf<fnref n="1" />, I. A. Grenier<fnref n="2" />, T. Joubaud<fnref n="2" />, and R. Schlickeiser<fnref n="1" />

Received 2019 April 15; revised 2019 May 21; accepted 2019 May 23; published 2019 July 8

## Abstract

A solution of the transport equation for cosmic rays in turbulent magnetic fields in a spherically symmetric geometry is presented. The results are applied to particle propagation in superbubbles. In the fully analytical calculation, various energy-loss processes are considered. From the distribution function of the cosmic-ray particles, the distribution for pions from continuous losses is computed. Folding with the appropriate cross section yields the gamma-ray distribution. It is shown that in the case of efficient acceleration the volume-integrated $\pi^0$-decay gamma-ray emission from the superbubble exhibits the characteristic hard differential number density spectrum $N(E_\gamma) \propto E_\gamma^{-0.29}$ at large gamma-ray energies, $E_\gamma \gg 70$ MeV. This prediction is useful for identifying efficient superbubble cosmic-ray accelerators from gamma-ray observations.

**Key words:** acceleration of particles – cosmic rays – gamma rays: ISM – ISM: bubbles

## 1. Introduction

Cosmic rays are a long-standing problem in astrophysics and to this day, not all of the processes involved in their acceleration and propagation are entirely understood. Particularly, the acceleration to relativistic energies in strong sources of non-thermal radiation and the evaluation of how the particles diffuse and affect the interstellar medium are of particular interest to modern astrophysics. Due to their influence on the thermodynamical and chemical states of the dense clouds that lead to star formation, cosmic rays are a major factor in the Galactic evolution (Wentzel 1971). For a long time, diffusive first-order Fermi processes at shock fronts have been regarded as a prime candidate for energy accumulation of cosmic rays. Observations with the Fermi Gamma-ray Space Telescope and Cerenkov telescopes have indeed identified strong emission regions near the shock fronts of supernova remnants. Relativistic shocks form during the interaction of relativistic and super-Alfvenic outflows with the ambient ionized interstellar medium (Gerbig & Schlickeiser 2011). Their presence in starburst regions adds additional pressure to the Galactic winds, which in turn spread into the halo of the galaxy (Breitschwerdt et al. 1991).

The gamma-ray observations of the Cygnus region with the Fermi-LAT have been used to analyze the total mass of the different gas phases and the cosmic-ray interactions with this mass (Ackermann et al. 2012). The analysis of this active region, which has recently formed hundreds of massive stars at a distance of 1.4 kpc, has uncovered a 50 pc wide cocoon of freshly accelerated cosmic rays (Ackermann et al. 2011), which flood the cavities created by the stellar winds and ionization fronts of the stellar clusters. These cosmic rays show a momentum spectrum that is significantly harder than the power-law spectrum that can be measured in this spiral arm after the propagation through the Galactic disk and halo. This recent discovery has disclosed an important new facet of the problem, What is the impact on particle (re-)acceleration and diffusion of the significant level of turbulence generated in the superbubble medium by the activity of the numerous massive stars?

Massive stars commonly cluster in space and time (Blaauw 1964). Giant O and B stars produce strong stellar winds that sweep the surrounding interstellar gas and push it outward. Winds and supernovae collectively power the turbulent ionized gas, magnetic fields, and shocks of the superbubble (Krause et al. 2014). The integrated energy released from a single wind is comparable with the energy yield of a supernova (Parizot et al. 2004; Cesarsky & Montmerle 1983).

For isolated stars, the deposited energy was examined closely (Seo et al. 2018) by decomposing the live cycle of the massive stars into different phases and using empirical data to model the radiation. It was found that the acceleration due to the stellar wind is significant, but the contribution from supernova explosions is still about four times larger. The processes that are involved in the transfer of mechanical energy of the wind into turbulence are numerous and difficult to describe accurately even for isolated systems. For a non-isolated system that includes colliding winds and subsequent supernova events, and therefore also collisions between supernova remnants and the individual stellar winds, the relative influence of these processes cannot be estimated quantitatively.

Due to their proximity, the stellar wind bubbles eventually merge and from their unification, a superbubble is formed. The supersonic winds and explosions generate fast shock waves and supersonic turbulence in the surrounding medium. Additional to the plasma instabilities generated by merging the stellar winds from the different OB stars and their associated magnetic and velocity fluctuations, the OB stars themselves can have strong magnetic fields $O(kG)$ (Castro et al. 2017).

Including all of the above mentioned influences, superbubbles appear as powerful sources of mechanically and magnetically turbulent flows where traversing particles are potentially scattered and accelerated. The individual and spatially separated sources of radiative pressure create a complex system of individual shock fronts and velocity fields. Additionally, supernova remnants that sweep the interior of the bubble region provide further acceleration capacity by picking up particles inside of the bubble. Thus, the particle diffusion and acceleration by repeated shocks may substantially modify the properties of cosmic rays traversing this region.

<footnote n="1">¹ Institut für Theoretische Physik, Lehrstuhl IV:Weltraum- und Astrophysik, Ruhr-Universität Bochum D-44780 Bochum, Germany; tt@tp4.rub.de, rsch@tp4.rub.de</footnote>

<footnote n="2">² Laboratoire AIM Paris-Saclay, CEA/Irfu, CNRS, Université Paris Diderot, F-91191 Gif-sur-Yvette Cedex, France; isabelle.grenier@cea.fr, theo.joubaud@cea.fr</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Model of the superbubble used in this paper: galactic cosmic rays penetrating the outer rim of the bubble region at $R_0$ are scattered and accelerated on the magnetic turbulence and the random velocity fields within the bubble. These fields confine the particles within the bubble and accelerate them stochastically. The physical extend of the bubble is denoted with $R_0$, from that distance onwards, the particles can stream away freely. For simplicity, the free streaming boundary condition is replaced using the free-escape boundary condition (19) instead. Particles that reach a certain distance $R_1$ from the center of the bubble take no longer part in the acceleration process.</figcaption></figure>

## 2. Particle Transport in Spherical Symmetry

In this paper, the superbubble is modeled as a spherical region with a high level of magnetic turbulence (Figure 1). The total cosmic-ray flux is determined by freshly accelerated cosmic rays generated in the center of the bubble and reaccelerated cosmic rays that are injected from the outside. These cosmic rays are assumed to be preaccelerated to a power law-like distribution function, similar to the distribution measured in the solar system.

The transport and acceleration of energetic particles in the partially turbulent magnetic fields associated with shocks is described using the Fokker-Planck equation for the particle distribution function (Schlickeiser 2011). The diffusion approximation for the particle density in the rest frame of the fluid is a well-known simplified form of the Fokker-Planck equation, which results when turbulent pitch-angle scattering is strong enough to ensure that the scale of the particle density variation is significantly greater than the particle mean free path (Jokipii 1966). Numerical studies confirmed the accuracy of the diffusion approximation in a uniform mean guide magnetic field (Kota et al. 1982).

### 2.1. Acceleration Processes in Superbubbles

In a medium, containing MHD plasma wave turbulence and moving with non-relativistic speed $U$, suprathermal and cosmic-ray particles obey the transport equation (Schlickeiser 2002, 2015)

$$ \frac{\partial f}{\partial t} + \mathcal{R}f + \nabla \cdot [Uf - K\nabla f] - \frac{1}{p^2} \frac{\partial}{\partial p} \left[ D_{pp,MHD} p^2 \frac{\partial f}{\partial p} + \frac{\nabla \cdot U}{3} p^3 f \right] = S(\mathbf{x}, p, t), \quad (1) $$

where

$$ \mathcal{R}f = p^{-2} \frac{\partial}{\partial p} [p^2 \dot{p}_{\text{loss}} f] + \frac{f}{T_c} \quad (2) $$

accounts for continuous ($\dot{p}_{\text{loss}}$) and catastrophic ($T_c$) momentum loss processes of cosmic-ray particles, and $S(\mathbf{x}, p, t)$ denotes the injection rate. $F$ denotes the omnidirectional particle distribution, $K$ is the spatial diffusion tensor, and $D_{pp} = a_1 p^2 / \kappa_{zz}$ is the momentum diffusion coefficient, which for MHD plasma wave turbulence is inversely proportional to the parallel spatial diffusion coefficient. Besides this stochastic acceleration, acceleration of particles arises from the term $\nabla \cdot U$. The latter term gives rise to first-order Fermi acceleration of charged particles in regions of space, where $\nabla \cdot U < 0$ (often referred to adiabatic compression), and to deceleration of particles in regions of space where $\nabla \cdot U > 0$ (referred to as adiabatic deceleration).

However, this term can also lead to second-order Fermi acceleration (described by an additional momentum diffusion coefficient) in physical systems containing random velocity fluctuations producing regions of space with both $\nabla \cdot U > 0$ and $\nabla \cdot U < 0$. This case has been investigated in detail before by Bykov & Toptygin (1982, 1993, 2001), Ptuskin (1988), Fedorov & Stehlik (2010), and Jokipii & Lee (2010). A suitable ensemble-averaging of the transport Equation (1) over the random velocity fluctuations in a frame of reference where $\langle U \rangle = 0$ then leads to (see Section 3 of Jokipii & Lee 2010 for mathematical details)

$$ \frac{\partial F}{\partial t} + \mathcal{R}F - \nabla \cdot K \nabla F - \frac{1}{p^2} \frac{\partial}{\partial p} p^2 \left[ (D_{pp,MHD} + D_{pp,VF}) \frac{\partial F}{\partial p} \right] = S(\mathbf{x}, p, t) \quad (3) $$

for the ensemble-averaged distribution function $F = \langle f \rangle$.

The transport Equation (3) also applies to the related problem of particle acceleration at multiple shocks separated by rarefactions, which has been considered before by Schneider (1993) and Melrose & Crouch (1997). Given their long ages, superbubbles harbor many stellar wind shocks with non-relativistic shock speeds, such that an ensemble-average over the acceleration process at different shocks is necessary. Clearly, the random velocity fluctuations associated with the shock statistics determine the additional momentum diffusion coefficient $D_{pp,VF}$ from velocity fluctuations in Equation (3).

### 2.2. Relevant Timescales

Let $F(r, p)$ be the cosmic-ray distribution function and $s(r, p)$ be the source function for the injection of cosmic rays for a spherical superbubble. The momentum losses are dominated by pion production and can be described as a continuous loss process with the rate (Schlickeiser 2009)

$$ \frac{\dot{p}_L}{p} = \frac{1}{\tau_L} = n_{\text{gas}} b_{\pi} s^{-1}, \quad (4) $$

where $b_{\pi} = 5.85 \times 10^{-16} \text{ cm}^3 \text{ s}^{-1}$ represents the pion-decay constant. With the typical high gas density $n_{\text{gas}} = 10n_1 \text{ cm}^{-3}$ inside the superbubble the momentum-independent pion
cooling time is

$$ \tau_L = \frac{1}{n_{\text{gas}} b_{\pi}} = 5.4 \cdot 10^6 n_1^{-1} \text{ yr.} \qquad (5) $$

Within this paper, we consider a spatially uniform diffusion coefficient

$$ k(r, p) = k_0 \chi(p) = 10^{25} k_{25} \chi(p) \text{ cm}^2 \text{ s}^{-1}. \qquad (6) $$

Its value is about an order of magnitude greater than the typical interstellar spatial diffusion coefficient accounting for the enhanced level of MHD turbulence inside the superbubble. Throughout this work we concentrate on the momentum diffusion $D_{pp,MHD}$ from particle-wave interactions with MHD plasma waves. The spatial diffusion coefficient and the momentum diffusion coefficient due to resonant interactions with MHD turbulence are closely related via (Schlickeiser 2002)

$$ \begin{gathered} a_1 = \frac{1}{\tau_A(p)} \approx \frac{v_A^2}{9 \chi(p) k_0} \approx 10^{-12} \frac{v_7^2}{k_{25} \chi(p)} \text{ s}^{-1}, \\ \tau_A(p) = 3.2 \cdot 10^4 \frac{k_{25} \chi(p)}{v_7^2} \text{ yr} \end{gathered} \qquad (7) $$

with the typical value of the Alfvén speed $v_A = 100v_7 \text{ km s}^{-1}$.

The third timescale of interest is the escape time of particles from the superbubble by spatial diffusion, which is given by

$$ \tau_D(p) = \frac{R_0^2}{k_0 \chi(p)} = 3.0 \cdot 10^4 \frac{R_0^2(\text{pc})}{k_{25} \chi(p)} \text{ yr}, \qquad (8) $$

which is of the same order as the acceleration timescale (7). The diffusion and acceleration timescales are about two orders of magnitude smaller than the pion loss time (5). It is convenient to introduce the ratio between acceleration and diffusion timescales

$$ \xi^2 = \frac{\tau_A}{\tau_D} = \left( \frac{k_{25} \chi(p)}{R_0 pc v_7} \right)^2, \qquad (9) $$

as a measure for the effectiveness of the accelerator.

Small values of $\xi$ correspond to an efficient acceleration of cosmic-ray particles, so the required time for the process is much smaller than the time spent inside the accelerator. Large values of $\xi$, however, correspond to an inefficient acceleration: the particles exit the bubble, before any efficient acceleration can occur. For prescribed injection conditions the steady-state cosmic-ray momentum spectrum is determined by the relative importance of these three timescales: (5), (7), and (8). The steady-state transport equation in the spherically symmetric superbubble, which accounts for momentum diffusion, pion losses, and spatial diffusion, can be written as

$$ \begin{gathered} \frac{1}{r^2} \frac{\partial}{\partial r} r^2 k_0 \frac{\partial}{\partial r} F + \frac{1}{p^2 \chi(p)} \frac{\partial}{\partial p} \left( \frac{a_1 p^4}{\chi(p)} \frac{\partial}{\partial p} + p^2 \dot{p}_L \right) F = - \frac{s(r, p)}{\chi(p)}. \end{gathered} \qquad (10) $$

$L_r$ represents the spatial diffusion of the particles inside the bubble, the first and second terms in $L_p$ account for stochastic

acceleration and the radiative losses, respectively, and both terms are balanced by the source rate $s(r, p)$.

## 2.3. Limiting Cases

The challenge is to determine the quasi-equilibrium momentum spectrum of cosmic rays resulting from the simultaneous action of momentum diffusion, pion momentum losses, and spatial diffusion on an arbitrary but given steady-state injection rate as described by Equation (10). Of course, it is straightforward (without any complicated mathematics) to directly infer the quasi-equilibrium momentum spectrum, if only one of these three processes is acting:

(i) If momentum diffusion is dominating, the equilibrium spectrum can be constructed from the condition that the diffusive flux in momentum space vanishes, i.e.,

$$ \frac{1}{p^2 \chi(p)} \frac{\partial}{\partial p} \left[ \frac{a_1 p^4}{\chi(p)} \frac{\partial F}{\partial p} \right] = 0, \qquad (11) $$

yielding for $\chi(p) = 1$ readily the two solutions $F_1(p) = \text{const.}$ and $F_2(p) = p^{-3}$. For the power-law dependence $\chi(p) \propto p^{\eta}$ one obtains $F_1(p) = \text{const.}$ and $F_2(p) = p^{\eta-3}$.

(ii) If pion momentum losses dominate one finds for the quasi-equilibrium momentum spectrum

$$ F(r, p) = \frac{\int_{p}^{\infty} dp' p'^2 s(r, p')}{p^2 \dot{p}_L}, \qquad (12) $$

such that for power-law injection $s(r, p) \propto p^{-k}$ the quasi-equilibrium momentum spectrum exhibits the same momentum dependence $\propto p^{-k}$.

(iii) For dominating spatial diffusion the quasi-equilibrium momentum spectrum is proportional to $s(r, p)/\chi(p)$.

However, such simple conclusions on the quasi-equilibrium momentum spectrum are no longer possible if all three processes operate simultaneously with comparable operation timescales as in our case. Then one has to employ a more involved mathematical apparatus, as we show next.

## 2.4. Solution of the Transport Equation

Using separation of variables, the spatial and momentum part of the equation can be decoupled. By defining

$$ F(r, p) = \sum_{m} f_{m}(p) H_{m}(r) \qquad (13) $$

and expanding the source term using the spatial eigenfunctions $H_m$

$$ s(r, p) = s_1(p)s_2(r) = s_1(p) \sum_{m} s_m H_m(r), \qquad (14) $$

the spatial part (Section 2.4.1) and the momentum part (Section 2.4.3) can be treated separately.

### 2.4.1. Spatial Distribution

The spatial eigenfunctions are given as the solutions of the homogeneous equation

$$ \frac{1}{r^2} \frac{d}{dr} r^2 k_0 \frac{d}{dr} H_m + \lambda_m H_m = 0. \qquad (15) $$
To further treat the spatial part, multiply Equation (15) with $r^4$ and substitute $r = 1/x$ providing

$$ \frac{d^2 H_m}{dx^2} + \frac{\lambda_m}{k_0} x^{-4} H_m = 0. \qquad (16) $$

This differential equation leads to Bessel functions with index (1/2) with the two solutions

$$ H_1 = \frac{\sin\left(\sqrt{\frac{\lambda_m}{k_0}} r\right)}{r}, \quad H_2 = \frac{\cos\left(\sqrt{\frac{\lambda_m}{k_0}} r\right)}{r}. \qquad (17) $$

Since the distribution function is supposed to be finite at the origin, only the first solution $H_1$ can contribute. It is possible to choose either the Marshak condition, which is equivalent to a free streaming of particles

$$ \left[ k_0 \frac{\partial F}{\partial r} - \frac{c}{2} F \right]_{r=R_0} = 0 \qquad (18) $$

at a distance where the superbubble turbulence is sufficiently damped. Alternatively, we use the simpler free-escape boundary condition (Lerche & Schlickeiser 1985),

$$ F(r = R_1, p) = 0, \qquad (19) $$

at

$$ R_1 = R_0 + \frac{2k_0\chi}{c} \approx R_0 + 0.002k_{25} \text{ pc}. \qquad (20) $$

With the free-escape boundary condition the eigenvalues $\lambda_m$ are fixed by the roots of the sine function

$$ H_m(r = R_1) \stackrel{!}{=} 0 \quad \Rightarrow \quad \lambda_m = \frac{m^2 \pi^2 k_0}{R_1^2}. \qquad (21) $$

with $m = 1, 2, 3, ...$. Substituting the eigenvalues into the spatial eigenfunction yields the spatial dependence of the distribution function inside the superbubble

$$ H_m(r) = r^{-\frac{1}{2}} J_2\left(m\pi \frac{r}{R_1}\right) \propto \frac{\sin\left(m\pi \frac{r}{R_1}\right)}{r}. \qquad (22) $$

### 2.4.2. Orthogonality Relations of the Spatial Functions

The eigenfunctions form a complete set of orthogonal functions. Multiply Equation (15) by $H_l$

$$ \frac{1}{r^2} \left( \frac{\partial}{\partial r} r^2 k_0 \frac{\partial}{\partial r} H_m \right) H_l + \lambda_m H_m H_l = 0, \qquad (23) $$

interchange the indices

$$ \frac{1}{r^2} \left( \frac{\partial}{\partial r} r^2 k_0 \frac{\partial}{\partial r} H_l \right) H_m + \lambda_l H_l H_m = 0, \qquad (24) $$

and subtract these equations. Integrating from 0 to $R_1$ yields the orthogonality relations

$$ 0 = \int_{0}^{R_1} dr \left\{ H_l \frac{\partial}{\partial r} r^2 k_0 \frac{\partial}{\partial r} H_m - H_m \frac{\partial}{\partial r} r^2 k_0 \frac{\partial}{\partial r} H_l + \frac{\lambda_m - \lambda_l}{k_0} r^2 H_m H_l \right\} \qquad (25) $$

$$ = \int_{0}^{R_1} dr \left\{ -r^2 \frac{\left( \frac{\partial H_l}{\partial r} \frac{\partial H_m}{\partial r} - \frac{\partial H_m}{\partial r} \frac{\partial H_l}{\partial r} \right)}{\underbrace{\frac{\partial H_l}{\partial r}}_{=0}} + \frac{\lambda_m - \lambda_l}{k_0} r^2 H_m H_l \right\} + \left[ H_l r^2 \frac{\partial H_m}{\partial r} - H_m r^2 \frac{\partial H_l}{\partial r} \right]_{0}^{R_1} \qquad (26) $$

$$ = \int_{0}^{R_1} dr \frac{\lambda_m - \lambda_l}{k_0} r^2 H_m H_l. \qquad (27) $$

Applying the eigenfunction decomposition (14), only the expansion coefficients remain to be calculated.

### 2.4.3. Momentum Distribution

Inserting the decomposition into spatial eigenfunctions to the original differential equation with the spatial and momentum operators $L_r$ and $L_p$, respectively, yields

$$ \begin{aligned} & L_r \sum_m f_m(p) H_m(r) + L_p \sum_m f_m(p) H_m(r) \\ & = - \sum_m s_m H_m(r) \frac{s_1(p)}{\chi(p)}. \end{aligned} \qquad (28) $$

The orthogonality relation is now used to simplify the differential equation. It acts on the sums, by reducing them to each element of the series individually. By integrating the equation with $\int_0^{R_1} dr r^2 H_l(r)$, all terms with $m \neq l$ vanish:

$$ \sum_m \left[ -f_m(p)\lambda_m + L_p f_m(p) + s_m \frac{s_1}{\chi(p)} \right] \int_0^{R_1} dr r^2 H_l H_m = 0. \qquad (29) $$

Thus,

$$ L_p f_m(p) - \lambda_m f_m(p) = -s_m \frac{s_1(p)}{\chi(p)}, \qquad (30) $$

such that only the momentum equation still remains to be solved:

$$ \frac{1}{p^2 \chi} \frac{d}{dp} \left( \frac{a_1 p^4}{\chi} \frac{d}{dp} + p^2 \dot{p}_L \right) f_m - \lambda_m f_m = -s_m \frac{s_1(p)}{\chi(p)}. \qquad (31) $$

## 2.5. Momentum-independent Spatial Diffusion Coefficient

To further advance the calculation, the spatial diffusion coefficient is assumed to be independent of momentum, i.e., $\chi(p) = 1$. With the pion loss rate (4) the remaining momentum equation reads

$$ \frac{1}{p^2} \frac{d}{dp} \left( p^4 \frac{d}{dp} + b p^3 \right) f_m - \frac{\lambda_m}{a_1} f_m = -s_m \frac{s_1(p)}{a_1} \qquad (32) $$

with<fnref n="3" />

$$ b = \frac{b_{\pi} n_{\text{gas}}}{a_1} = \frac{\tau_{\text{A}}}{\tau_{\text{L}}} = 5.9 \cdot 10^{-3} \frac{k_{25} n_1}{v_7^2}. \qquad (33) $$

<footnote n="3">³ The smallness of $b$ justifies ignoring the influence of pion losses completely by setting $b = 0$ in the transport Equation (32). We do not follow this approximation here, as Equation (32) can be analytically solved for finite values of $b$.</footnote>
The analytical solution to this equation is detailed in the Appendix. In terms of the proton Lorentz factor $p = m_p c \sqrt{\gamma^2 - 1}$ the differential particle density distribution reads

$$ 
\begin{aligned}
N(r, \gamma) &= \frac{m_p c \gamma}{\sqrt{\gamma^2 - 1}} N(r, p = m_p c \sqrt{\gamma^2 - 1}) \\
&= \frac{2\pi(m_p c)^{\frac{3-b}{2}}}{a_1 r} \gamma(\gamma^2 - 1)^{-\frac{1+b}{4}} \\
&\quad \times \sum_{m=1}^{\infty} \frac{s_m \sin\left(m\pi \frac{r}{R_l}\right)}{\mu_m} \left[ p^{-\mu_m} \int_0^p dp_0 s_1(p_0) p_0^{\frac{b+1}{2}+\mu_m} \right. \\
&\qquad \left. + p^{\mu_m} \int_p^{\infty} dp_0 s_1(p_0) p_0^{\frac{b+1}{2}-\mu_m} \right]_{p=m_p c \sqrt{\gamma^2-1}}.
\end{aligned} 
\quad (34) $$

## 2.6. Impact of the Acceleration Efficiency

The ratio of the acceleration and diffusion timescale plays a significant role in the resulting spectral index of the accelerated cosmic-ray distribution. In the following, the extreme cases of very efficient and very inefficient acceleration are discussed. Assuming a power-law spectrum with a strict lower limit $p_i = 10 m_p c$ as the momentum dependence of the injection distribution,

$$ s_1(p_0) = q_1 p_0^{-s} \theta(p_0 - p_i), \quad (35) $$

where a typical value for the spectral index of $s = 4.7$ is used. We compute the integral

$$ 
\begin{aligned}
& p^{-\mu_m} \int_0^p dp_0 s_1(p_0) p_0^{\frac{b+1}{2}+\mu_m} + p^{\mu_m} \int_p^{\infty} dp_0 s_1(p_0) p_0^{\frac{b+1}{2}-\mu_m} \\
&= q_1 \left\{ \begin{aligned} & \frac{p_i^{\frac{b+3}{2}-s}}{\mu_m + s - \frac{b+3}{2}} \left( \frac{p}{p_i} \right)^{\mu_m} && \text{for } p < p_i \\ & p^{\frac{b+3}{2}-s} \left[ \frac{1}{\mu_m + s - \frac{b+3}{2}} + \frac{1 - \left( \frac{p}{p_i} \right)^{s-\mu_m - \frac{b+3}{2}}}{\mu_m + \frac{b+3}{2} - s} \right] && \text{for } p \ge p_i. \end{aligned} \right.
\end{aligned} 
\quad (36) $$

Thus, the differential particle density distribution (14) can be expressed as

$$ 
\begin{aligned}
N(r, p) = & \frac{2\pi q_1 p^{\frac{1-b}{2}}}{a_1 r} \sum_{m=1}^{\infty} \frac{s_m \sin\left(m\pi \frac{r}{R_l}\right)}{\mu_m} \\
& \times \left\{ \begin{aligned} & \frac{p_i^{\frac{b+3}{2}-s}}{\mu_m + s - \frac{b+3}{2}} \left( \frac{p}{p_i} \right)^{\mu_m} && \text{for } p < p_i \\ & p^{\frac{b+3}{2}-s} \left[ \frac{1}{\mu_m + s - \frac{b+3}{2}} + \frac{1 - \left( \frac{p}{p_i} \right)^{s-\mu_m - \frac{b+3}{2}}}{\mu_m + \frac{b+3}{2} - s} \right] && \text{for } p \ge p_i. \end{aligned} \right.
\end{aligned} 
\quad (37) $$

Rewriting the momentum dependence in terms of the ratio between the injection threshold and the momentum $x = \frac{p}{p_i}$ and defining $\eta_m$ as shorthand for the $x$-dependent part of the

distribution function yields

$$ 
\begin{aligned}
N(r, p) = & \frac{2\pi q_1 p_i^{2-s}}{a_1 r} \sum_{m=1}^{\infty} \frac{s_m \sin\left(m\pi \frac{r}{R_l}\right)}{\mu_m} \eta_n(x), \quad \eta_m(x) \\
& = \begin{cases} \frac{x^{\mu_m + \frac{1-b}{2}}}{\mu_m + s - \frac{b+3}{2}} & \text{for } x < 1 \\ x^{2-s} \left[ \frac{1}{\mu_m + s - \frac{b+3}{2}} + \frac{1 - x^{s-\mu_m - \frac{b+3}{2}}}{\mu_m + \frac{b+3}{2} - s} \right] & \text{for } x \ge p_i \end{cases} 
\end{aligned} 
\quad (38) $$

with $\mu_m$ given in Equation (73). For inefficient accelerators, i.e., large values of $\xi \gg 1$, the parameter $\mu_m$ can be approximated as

$$ \mu_m \stackrel{(\xi \gg 1)}{=} \sqrt{\left(\frac{3-b}{2}\right)^2 + m^2 \xi^2} \approx m\xi \gg 1. \quad (39) $$

In the domain of $x < 1$ of Equation (38), the exponent is dominated by $m\xi$ so that the other terms in the exponent can be neglected. This results into a steeply rising spectrum, where almost no particles are present below the injection threshold. The second case $x \ge 1$ simplifies, as the term proportional to $x^{-\xi m}$ vanishes. The resulting distribution resembles the injection distribution. This can be understood as the particles leave the acceleration region before any acceleration could take place. The dependence of $\eta(x)$ for large $\xi$ then reduces to

$$ \eta_m(x) \approx \frac{1}{m\xi} \begin{cases} x^{\xi m} & \text{for } x < 1 \\ 2x^{2-s} & \text{for } x \gg 1. \end{cases} \quad (40) $$

For very efficient accelerators, the diffusive time is much larger than the acceleration time, which corresponds to small values of $\xi$. Therefore, the parameter $\mu_m$ is mostly determined by the first term

$$ \mu_m \stackrel{(\xi \ll 1)}{=} \sqrt{\left(\frac{3-b}{2}\right)^2 + m^2 \xi^2} \approx \frac{3-b}{2}. \quad (41) $$

In this limiting case, Equation (38) reads

$$ \eta_m(x) \approx \begin{cases} \frac{x^{2-b}}{s-b} & \text{for } x < 1 \\ x^{2-s} \left( \frac{1}{s-b} + \frac{1-x^{s-3}}{3-s} \right) & \text{for } x \ge 1. \end{cases} \quad (42) $$

For $x \gg 1$ and $s > 3$, the last term dominates and $\eta_m(x)$ reduces to

$$ \eta_m(x) \approx \frac{x^{-1}}{s-3}. \quad (43) $$

Efficient accelerators can result in a power spectrum as hard as $x^{-1}$ for large momenta and $x^2$ for energies below the injection threshold (Figure 2), in excellent agreement with the discussion in Section 2.3.

## 3. Pion Distribution

The neutral-pion production function from inelastic hadron-hadron collisions of cosmic rays in the superbubble can be calculated from the cosmic-ray density distribution (34) convolving $N(r, \gamma)$ with the pion production power of a
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Distribution function plotted against the momentum of the cosmic-ray particles. The blue curve depicts the distribution after acceleration in an ideal accelerator (small $\xi$). Furthermore, the distribution of the particles for the opposite case, i.e., an inefficient accelerator (large $\xi$), is represented by the orange line. The reaccelerated spectrum in the latter case is almost identical to the injection spectrum, while it is significantly harder for an efficient accelerator.</figcaption></figure>

single proton (Schlickeiser 2002)

$$ 
\begin{aligned}
P_{\pi}(\gamma, \gamma_{\pi}) &= 1.3c\gamma_{\pi}n_{\text{gas}}\sum_{pp}^{\pi_{\pi}}(\gamma)\delta(\gamma_{\pi} - \gamma_{\pi}^{3})\Theta(\gamma - \gamma_{\text{thr}}) \\
&= \frac{4 \cdot 1.3}{3}cn_{\text{gas}}\gamma_{\pi}\gamma_{\pi}^{\frac{1}{4}}\sum_{pp}^{\pi_{\pi}}(\gamma)\delta(\gamma - \gamma_{\pi}^{\frac{4}{3}})\Theta(\gamma - \gamma_{\text{thr}})
\end{aligned}
\quad (44) $$

with $\gamma_{\text{thr}} = 1.43$ and (Stecker 1973)

$$ 
\begin{aligned}
\sum_{pp}^{\pi_{\pi}}(\gamma) &= \xi\sigma_{pp}^{\pi_{\pi}}(\gamma) \\
&= \begin{cases} 6.13 \cdot 10^{-26}(\gamma - 1)^{7.64} \text{ cm}^2 & \text{for } 1.43 \le \gamma \le 1.75 \\ 8.12 \cdot 10^{-27}(\gamma - 1)^{0.53} \text{ cm}^2 & \text{for } 1.75 < \gamma. \end{cases}
\end{aligned}
\quad (45) $$

The neutral-pion production spectrum then reads

$$ Q_{\pi}(r, \gamma_{\pi}) = \frac{1.26}{\gamma_{\pi}m_{\pi}c^2} \int_{1}^{\infty} d\gamma \, P_{\pi}(\gamma, \gamma_{\pi}) N(r, \gamma) \quad (46) $$

$$ = \frac{4 \cdot 1.64}{3} \frac{cn_{\text{gas}}}{m_{\pi}c^2} \int_{1.43}^{\infty} d\gamma \, \gamma^{\frac{1}{4}} \sum_{pp}^{\pi_{\pi}^{0}}(\gamma) N(r, \gamma) \delta(\gamma - \gamma_{\pi}^{\frac{4}{3}}) \quad (47) $$

$$ = \frac{2.19cn_{\text{gas}}}{m_{\pi}c^2} \gamma_{\pi}^{\frac{1}{3}} N(r, \gamma_{\pi}^{\frac{4}{3}}) \sum_{pp}^{\pi_{\pi}^{0}} (\gamma_{\pi}^{\frac{4}{3}}) \Theta(\gamma_{\pi} - 1.31). \quad (48) $$

Using Equation (34) we obtain

$$ 
\begin{aligned}
Q_{\pi}(r, \gamma_{\pi}) = & \frac{4.38\pi m_{\pi} n_{\text{gas}}}{a_{1} r m_{\pi}} (m_{p} c)^{\frac{1-b}{2}} \gamma_{\pi}^{\frac{5}{3}} (\gamma_{\pi}^{\frac{8}{3}} - 1)^{-\frac{1+b}{4}} \\
& \times \sum_{pp}^{\pi_{\pi}^{0}} (\gamma_{\pi}^{\frac{4}{3}}) \Theta(\gamma_{\pi} - 1.31) \sum_{m=1}^{\infty} \frac{s_{m} \sin\left(m\pi \frac{r}{R_{1}}\right)}{\mu_{m}} \\
& \times \left[ p^{-\mu_{m}} \int_{0}^{p} dp_{0} s_{1}(p_{0}) p_{0}^{\frac{b+1}{2}+\mu_{m}} \right. \\
& \left. + p^{\mu_{m}} \int_{p}^{\infty} dp_{0} s_{1}(p_{0}) p_{0}^{\frac{b+1}{2}-\mu_{m}} \right]_{p=m_{p}c \sqrt{\frac{8}{\gamma_{\pi}^{2}}-1}}.
\end{aligned}
\quad (49) $$

In order to compare with volume-integrated $\gamma$-ray spectra from superbubbles, we need to calculate the volume-integrated

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Proton distribution (13) in equilibrium for galactic cosmic rays injected at the edge of the superbubble with the distribution function (35). The spectral index of the source function phase space density $s = 4.7$ was chosen to resemble the observed cosmic-ray distribution. The plotted function illustrates the integrated distribution function with energies between 10 GeV and 10 TeV. The normalization was chosen such that the density is set to be close to the boundary of the bubble. Both cases, efficient and inefficient acceleration, are presented here. The efficient accelerator is shown in blue with $\xi = 0.1$, while the inefficient accelerator is depicted as orange with $\xi = 10$. In the case of efficient acceleration, the cosmic-ray density first increases toward the center, because the particles close to the boundary can easily leave the bubble, until it reaches a maximum. From there on, the effect of the enhanced turbulence prevails, making it difficult for the particles to penetrate further inside the bubble. The impact of inefficient accelerators on the spatial distribution is much smaller. Since the level of turbulence is small, the particles can easily penetrate into the deeper layers of the bubble. The density at the edge is smaller compared to the center, as the particles close to the edge can stream away.</figcaption></figure>

neutral-pion production spectrum

$$ Q_{\pi}(\gamma_{\pi}) = 4\pi \int_{0}^{R_{0}} dr \, r^{2} Q_{\pi}(r, \gamma_{\pi}). \quad (50) $$

With the integral

$$ 
\begin{aligned}
& \int_{0}^{R_{0}} dr \, r \sin\left(m\pi \frac{r}{R_{1}}\right) \\
= & \frac{R_{1}^{2}}{m^{2}\pi^{2}} \left[ \sin\left(m\pi \frac{R_{0}}{R_{1}}\right) - m\pi \frac{R_{0}}{R_{1}} \cos\left(m\pi \frac{R_{0}}{R_{1}}\right) \right].
\end{aligned}
\quad (51) $$

we find

$$ 
\begin{aligned}
Q_{\pi}(\gamma_{\pi}) = & \frac{17.52m_{p}n_{\text{gas}}R_{1}^{2}(m_{p}c)^{\frac{1-b}{2}}}{a_{1}m_{\pi}} \gamma_{\pi}^{\frac{5}{3}}(\gamma_{\pi}^{\frac{8}{3}} - 1)^{-\frac{1+b}{4}} \\
& \times \sum_{pp}^{\pi_{\pi}^{0}} (\gamma_{\pi}^{\frac{4}{3}}) \Theta(\gamma_{\pi} - 1.31) \sum_{m=1}^{\infty} \frac{s_{m}}{m^{2}\mu_{m}} \left[ \sin\left(m\pi \frac{R_{0}}{R_{1}}\right) \right. \\
& \left. - m\pi \frac{R_{0}}{R_{1}} \cos\left(m\pi \frac{R_{0}}{R_{1}}\right) \right] \left[ p^{-\mu_{m}} \int_{0}^{p} dp_{0} s_{1}(p_{0}) p_{0}^{\frac{b+1}{2}+\mu_{m}} \right. \\
& \left. + p^{\mu_{m}} \int_{p}^{\infty} dp_{0} s_{1}(p_{0}) p_{0}^{\frac{b+1}{2}-\mu_{m}} \right]_{p=m_{p}c \sqrt{\frac{8}{\gamma_{\pi}^{2}}-1}}.
\end{aligned}
\quad (52) $$

## 4. Reacceleration of Galactic Cosmic Rays

In this section, the effect of turbulent stochastic reacceleration inside the superbubble on a power-law spectrum of galactic cosmic rays passing by and penetrating into the bubble is studied.
Considering galactic cosmic rays as a primary injection source, the spatial part of the source function can be represented as a delta distribution at the boundary of the bubble:

$$s_2(r) = \delta(r - R_0) = \sum_l s_l H_l(r). \qquad (53)$$

Using the orthogonality relation (27) on the source term, the series coefficients can be computed as

$$\int_{0}^{R_1} dr \, r^2 H_m(r) \delta(r - R_0) = \int_{0}^{R_1} dr \, r^2 H_m(r) \sum_{l} s_l H_l(r), \quad (54)$$

We use the injection spectrum discussed in Section 2.6, which was introduced as a power-law spectrum with a strict lower limit for the source function

$$s_1(p_0) = q_1 p_0^{-s} \theta(p_0 - p_i). \qquad (58)$$

The relevant integrals for the convolution with Green's function are already computed in that section and consequently the neutral-pion production spectrum (57) becomes

$$Q_{\pi}(\gamma_{\pi}) = \frac{24.78 \pi m_p n_{\text{gas}} q_1 R_0 R_1^{\frac{1}{2}} (m_p c)^{\frac{1-b}{2}}}{a_1 m_{\pi}} \gamma_{\pi}^{\frac{5}{3}} (\gamma_{\pi}^{\frac{8}{3}} - 1)^{-\frac{1+b}{4}} \sum_{pp} \pi_{pp}^{0} (\gamma_{\pi}^{\frac{4}{3}}) \Theta(\gamma_{\pi} - 1.31) p_{l}^{\frac{b+3}{2}-s} \sum_{m=1}^{\infty} \frac{\sin\left(m\pi \frac{R_0}{R_l}\right)}{m^2 \mu_m} \\ \times \left[ \sin\left(m\pi \frac{R_0}{R_l}\right) - m\pi \frac{R_0}{R_l} \cos\left(m\pi \frac{R_0}{R_l}\right) \right] \left\{ \frac{1}{\mu_m + s - \frac{b+3}{2}} \left( \frac{\gamma_{\pi}^{\frac{8}{3}} - 1}{\gamma_{\pi,i}^{\frac{8}{3}} - 1} \right)^{\frac{\mu_m}{2}} \right. \\ \left. \qquad \times \left[ \frac{1}{\mu_m + s - \frac{b+3}{2}} + \frac{1 - \left( \frac{\gamma_{\pi}^{\frac{8}{3}} - 1}{\gamma_{\pi,i}^{\frac{8}{3}} - 1} \right)^{\frac{2s-2\mu_m-b-3}{4}}}{\mu_m + \frac{b+3}{2} - s} \right] \right\} \quad \begin{array}{l} \text{for } \gamma_{\pi} \le \gamma_{\pi,i} \\ \\ \text{for } \gamma_{\pi} \ge \gamma_{\pi,i} \end{array} \qquad (59)$$

yielding for the spatial expansion coefficient

with

$$s_m = \frac{R_0^2 H_m(R_0)}{\int_0^{R_1} dr r^2 H_m^2(r)} = \frac{2 R_0^{\frac{3}{2}} J_1\left(m \pi \frac{R_0}{R_1}\right)}{R_1^2 J_2^{\frac{3}{2}}(m\pi)} \\ = \frac{\pi \sqrt{2mR_0}}{R_1^{\frac{3}{2}}} \sin\left(m\pi \frac{R_0}{R_1}\right), \qquad (55)$$

$$\gamma_{\pi,i} = \left[ 1 + \left( \frac{p_i}{m_p c} \right)^2 \right]^{\frac{3}{2}}. \qquad (60)$$

where we used

$$\int_{0}^{R_1} dr r^2 H_m^2(r) = \int_{0}^{R_1} dr \, r J_2^2\left(m\pi \frac{r}{R_1}\right) \\ = \frac{1}{2} R_1^2 J_2^2(m\pi). \qquad (56)$$

With this value of $s_m$ the neutral-pion production spectrum (52) becomes

$$Q_{\pi}(\gamma_{\pi}) = \frac{24.78 \pi m_p n_{\text{gas}} R_0 R_1^{\frac{1}{2}} (m_p c)^{\frac{1-b}{2}}}{a_1 m_{\pi}} \gamma_{\pi}^{\frac{5}{3}} (\gamma_{\pi}^{\frac{8}{3}} - 1)^{-\frac{1+b}{4}} \\ \times \sum_{pp} \pi_{pp}^{0} (\gamma_{\pi}^{\frac{4}{3}}) \Theta(\gamma_{\pi} - 1.31) \sum_{m=1}^{\infty} \frac{\sin\left(m\pi \frac{R_0}{R_1}\right)}{m^2 \mu_m} \left[ \sin\left(m\pi \frac{R_0}{R_1}\right) \right. \\ \left. - m\pi \frac{R_0}{R_1} \cos\left(m\pi \frac{R_0}{R_1}\right) \right] \left[ p^{-\mu_m} \int_0^p dp_0 s_1(p_0) p_0^{\frac{b+1}{2}+\mu_m} \right. \\ \left. + p^{\mu_m} \int_p^{\infty} dp_0 s_1(p_0) p_0^{\frac{b+1}{2}-\mu_m} \right]_{p=m_p c \sqrt{\gamma_{\pi}^{\frac{8}{3}}-1}} \qquad (57)$$

The steady state spatial distribution function is shown in Figure 3 for the case of efficient and inefficient acceleration. Figure 4 depicts the spatially integrated distribution as a function of momentum for the cosmic-ray protons and the neutral pions.

### 4.1. Gamma-Ray spectrum

Neutral pions decay almost exclusively into two photons, therefore the γ-ray distribution function is closely related to the pion distribution (Schlickeiser 2002)

$$Q_{\gamma}(E_{\gamma}) = 2 \int_{\frac{E_{\gamma}}{m_{\pi}c^2} + \frac{m_{\pi}c^2}{4E_{\gamma}}}^{\infty} d\gamma_{\pi} \frac{Q_{\pi}(\gamma_{\pi})}{\sqrt{\gamma_{\pi}^2 - 1}}. \qquad (61)$$

Equation (61) can be used to determine the universal spectral index for the gamma-ray distribution. For large values of $\gamma_{\pi}$, the pion distribution function becomes a power law

$$Q_{\pi}(\gamma_{\pi}) \propto \gamma_{\pi}^{-s_{\pi}}. \qquad (62)$$

Energies far below the pion rest energy $E_{\gamma} \ll m_{\pi}c^2$ and far above $E_{\gamma} \gg m_{\pi}c^2$ consequently result in a gamma-ray spectrum that resembles a power law with the same spectral index as the pion distribution. This can be seen from the lower boundary of the integration, which becomes large at these
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left panel: the gamma-ray distribution from pion decay as a function of energy. Both cases (efficient and inefficient accelerator) are normalized to the value of the distribution function at the pion mass. Additionally, the distribution function for the inefficient accelerator is divided by four to provide better distinguishability. The pion distribution function below the injection energy is orders of magnitude smaller than that above $E_i$ and thus the contribution to the gamma-rays is negligible. This leads to the formation of the plateau in the diagram. Beyond the plateau region, the distribution transitions into a power-law shape. The spectral index in that regime coincides with the spectral index of the pion distribution. Right panel: the spectral energy distribution of the gamma-rays from pion decay. It has the same normalization as the first diagram. The efficient and inefficient accelerators separate especially well at high energies. While the total energy emitted into high-energy photons increases with higher energies for efficient accelerators, an inefficient accelerator cannot provide sufficient amounts of energy, which in turn results in less total energy radiated in high-energy photons.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Spatially integrated momentum distribution of the protons (left side) and pions (right side) as a function of the Lorenz factor inside the superbubble. The color coding from Figure 3 is consistently used for these diagrams as well. The high-energy behavior resembles a power law. In the low-energy region, the dominant feature of the function is a break at the injection momentum threshold. In the pion distribution, the shape of the production cross section (45) factors in which results in an additional break due to the lower cross section in that energy domain. This flattening of the distribution is only visible for the efficient accelerator case, as opposed to the inefficient case where the distribution function itself already rises steeply and thus the break is below the displayed area.</figcaption></figure>

energies

$$ \frac{E_{\gamma}}{m_{\pi}c^2} + \frac{m_{\pi}c^2}{4E_{\gamma}} \gg 1. \qquad (63) $$

Therefore, the one in the denominator can be dropped and the integral reduces to

$$ \begin{aligned} Q_{\gamma}(E_{\gamma}) &\propto 2 \int_{\frac{m_{\pi}c^2}{E_{\gamma}} + \frac{m_{\pi}c^2}{4E_{\gamma}}}^{\infty} d\gamma_{\pi} \frac{\gamma_{\pi}^{-s_{\pi}}}{\gamma_{\pi}} = 2 \int_{\frac{m_{\pi}c^2}{E_{\gamma}} + \frac{m_{\pi}c^2}{4E_{\gamma}}}^{\infty} d\gamma_{\pi} \gamma_{\pi}^{-s_{\pi}-1} \\ &= \frac{2}{s_{\pi}} \begin{cases} \left( \frac{4E_{\gamma}}{m_{\pi}c^2} \right)^{s_{\pi}} & \text{for } E_{\gamma} \ll m_{\pi}c^2 \\ \left( \frac{E_{\gamma}}{m_{\pi}c^2} \right)^{-s_{\pi}} & \text{for } E_{\gamma} \gg m_{\pi}c^2. \end{cases} \end{aligned} \qquad (64) $$

The hardest gamma-ray spectrum that can be produced by an efficient accelerator with $\xi \ll 1$ will therefore follow a power-law distribution in the tail regions with a spectral index of

$$ s_{\gamma} \stackrel{\xi \to 0}{=} \frac{4}{3} - 1.04 \approx 0.29. \qquad (65) $$

In the case of an inefficient accelerator, the high-energy gamma-ray spectral index approaches

$$ s_{\gamma} \stackrel{\xi \to \infty}{=} \frac{4}{3} \cdot 2.7 - 1.04 \approx 2.56, \qquad (66) $$

which is significantly greater than 0.29 in the case of an efficient accelerator. The gamma-ray distributions for both cases, efficient and inefficient acceleration, are shown in Figure 5. Particularly, the high-energy behavior of the distribution function serves a probe for the acceleration efficiency $\xi$.

## 5. Summary and Conclusions

Superbubbles are a potential source for efficient cosmic-ray acceleration due to the enhanced level of MHD turbulence reaccelerating penetrating cosmic rays from the ambient interstellar medium. The importance of this mechanism can be investigated by calculating the flux of gamma-rays from the decay of neutral pions generated in inelastic hadron–hadron collisions of the reaccelerated cosmic-ray hadrons with the superbubble gas. In this work we calculate the volume-integrated
gamma-ray spectrum from a spherical superbubble, accounting for the competing effects of pion cooling, stochastic reacceleration and particle escape by spatial diffusion. The steady-state cosmic-ray transport equation is solved analytically for the case in which interstellar cosmic rays with a power-law momentum distribution function enter the superbubble. The momentum spectrum of reaccelerated particles differs from the injection distribution, implying spectral differences of the pion-decay gamma-ray energy spectrum from superbubbles as compared to the gamma-ray energy spectrum from the outside interstellar medium.

The results of the analytic calculation are applied to the edge cases of potential accelerators. The momentum spectrum of reaccelerated cosmic rays sensitively depends on the parameter $\xi = \sqrt{\frac{\tau_A}{\tau_D}}$ determined by the ratio between the acceleration and diffusion timescales. Accelerators, which are unable to contain the particles for a sufficiently large time to amount in acceleration, only result in inefficient acceleration $\xi \ll 1$ for the cosmic rays and thus no spectral differences to the injection spectrum are found. Alternatively, for the case of efficient acceleration $\xi \ll 1$ approaches the behavior $N(p > p_i) \propto p^{-1}$, whereas for $N(p \le p_i) \propto p^{2-b}$. This also implies significant spectral differences for the pion decay and gamma-ray energy spectra as compared to the gamma-ray energy spectra from the ambient interstellar medium. The universal spectral index for the gamma-ray distribution was found to be $s_\gamma = 0.29$ for the case of extremely efficient acceleration. These predictions can be used to identify superbubbles with the required properties to efficiently accelerate cosmic rays. Detecting the implied differences in the volume-integrated gamma-ray spectrum with the observed gamma-ray spectra from superbubbles indicates that superbubbles are an important source of efficiently accelerating cosmic rays.

This work was financially supported by the collaborative France-German ANR-DFG project (Schl 201/34-1 and ANR-15-CE31-0019-01).

# Appendix
## Derivation of the Momentum Solution to the Transport Equation

In this appendix, we detail the derivation of the solution to the momentum part of the differential equation for arbitrary source distributions $s(r, p)$ by calculating the Green's function of the equation. The goal of this section is to derive an analytical solution to the remaining momentum part of the transport Equation (32). Using the ansatz

$$f_m(p) = p^{-b}M(p) \qquad (67)$$

provides for $M(p)$ the self-adjoint equation

$$\frac{d}{dp} \left( p^{4-b} \frac{dM}{dp} \right) - \frac{\lambda_m}{a_1} p^{2-b} M(p) = - \frac{s_m}{a_1} p^2 s_1(p), \quad (68)$$

which is solved for arbitrary injection rates $s_1(p)$ as

$$M(p) = \frac{s_m}{a_1} \int_0^\infty dp_0 G(p, p_0) p_0^2 s_1(p_0) \quad (69)$$

in terms of Green's function satisfying

$$\frac{d}{dp} \left( p^{4-b} \frac{dG}{dp} \right) - \frac{\lambda_m}{a_1} p^{2-b} G = -\delta(p - p_0). \quad (70)$$

### A.1. Green's Function

Green's function is computed from the solutions of the homogeneous equation corresponding to Equation (70)

$$h_1(p) = p^{\alpha_1}, \quad h_2(p) = p^{\alpha_2} \qquad (71)$$

with

$$\alpha_1 = \mu_m - \frac{3-b}{2}, \quad \alpha_2 = -\mu_m - \frac{3-b}{2}, \\ \mu_m = \sqrt{\left(\frac{3-b}{2}\right)^2 + \frac{\lambda_m}{a_1}}. \qquad (72)$$

We note that

$$\mu_m = \sqrt{\left(\frac{3-b}{2}\right)^2 + \frac{m^2 \pi^2 k_0}{R_1^2 a_1}} = \sqrt{\left(\frac{3-b}{2}\right)^2 + m^2 \xi^2}, \quad (73)$$

using (9). Obviously, $h_1(p)$ vanishes for $p \to 0$, whereas $h_1(p)$ vanishes for $p \to \infty$. We construct Green's function as

$$G(p, p_0) = \frac{1}{A} \begin{cases} h_1(p) h_2(p_0) & \text{for } p \le p_0 \\ h_1(p_0) h_2(p) & \text{for } p \ge p_0 \end{cases} \\ = \frac{1}{A} \begin{cases} p^{\alpha_1} p_0^{\alpha_2} & \text{for } p \le p_0 \\ p_0^{\alpha_1} p^{\alpha_2} & \text{for } p \ge p_0. \end{cases} \qquad (74)$$

The normalization factor $A$ is obtained by integrating the differential Equation (70) for the Green's function

$$\lim_{\eta \to 0} \int_{p_0 - \eta}^{p_0 + \eta} dp \times \left[ \frac{d}{dp} p^{4-b} \frac{dG}{dp} - \frac{\lambda_m}{a_1} p^{2-b} G + \delta(p - p_0) \right] = 0. \quad (75)$$

Since $G$ is continuous by construction, the second term vanishes. Inserting Equation (74) yields

$$A = -\lim_{\eta \to 0} \left[ p^{4-b} \frac{dG}{dp} \bigg|_{p_0+\eta} - p^{4-b} \frac{dG}{dp} \bigg|_{p_0-\eta} \right] \\ = (\alpha_1 - \alpha_2) p_0^{3-b+\alpha_1+\alpha_2} = 2\mu_m. \qquad (76)$$

Consequently we obtain

$$G(p, p_0) = \frac{(pp_0)^{\frac{b-3}{2}}}{2\mu_m} \begin{cases} (p/p_0)^{\mu_m} & \text{for } p \le p_0 \\ (p/p_0)^{-\mu_m} & \text{for } p \ge p_0. \end{cases} \qquad (77)$$
## A.2. Solution to the Momentum Problem

With Equations (67) and (69) we obtain the solution

$$f_m(p) = \frac{s_m}{2a_1\mu_m} p^{-\frac{b+3}{2}} \left[ p^{-\mu_m} \int_0^p dp_0 s_1(p_0) p_0^{\frac{b+1}{2}+\mu_m} + p^{\mu_m} \int_p^\infty dp_0 s_1(p_0) p_0^{\frac{b+1}{2}-\mu_m} \right], \quad (78)$$

implying, according to Equation (13),

$$F(r, p) = \frac{p^{-\frac{b+3}{2}}}{2a_1r} \sum_{m=1}^{\infty} \frac{s_m \sin\left(m\pi \frac{r}{R_1}\right)}{\mu_m} \left[ p^{-\mu_m} \times \int_0^p dp_0 s_1(p_0) p_0^{\frac{b+1}{2}+\mu_m} + p^{\mu_m} \int_p^\infty dp_0 s_1(p_0) p_0^{\frac{b+1}{2}-\mu_m} \right]. \quad (79)$$

The differential particle density distribution is then given by

$$N(r, p) = 4\pi p^2 F(r, p) = \frac{2\pi p^{\frac{1-b}{2}}}{a_1 r} \sum_{m=1}^{\infty} \frac{s_m \sin\left(m\pi \frac{r}{R_1}\right)}{\mu_m} \times \left[ p^{-\mu_m} \int_0^p dp_0 s_1(p_0) p_0^{\frac{b+1}{2}+\mu_m} + p^{\mu_m} \int_p^\infty dp_0 s_1(p_0) p_0^{\frac{b+1}{2}-\mu_m} \right]. \quad (80)$$

## ORCID iDs

T. Tolksdorf https://orcid.org/0000-0001-7506-7275

R. Schlickeiser https://orcid.org/0000-0003-3171-5079

## References

* Ackermann, M., Ajello, M., Allafort, A., et al. 2011, *Sci*, 334, 1103

* Ackermann, M., Ajello, M., Allafort, A., et al. 2012, *A&A*, 538, A71

* Blaauw, A. 1964, *ARA&A*, 2, 213

* Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, *A&A*, 245, 79

* Bykov, A. M., & Toptygin, I. N. 1982, *JGZG*, 50, 221

* Bykov, A. M., & Toptygin, I. N. 1993, *PhyU*, 36, 1020

* Bykov, A. M., & Toptygin, I. N. 2001, *AstL*, 27, 625

* Castro, N., Fossati, L., Hubrig, S., et al. 2017, *A&A*, 597, L6

* Cesarsky, C. J., & Montmerle, T. 1983, *SSRv*, 36, 173

* Fedorov, Y., & Stehlik, M. 2010, *JPhB*, 43, 185701

* Gerbig, D., & Schlickeiser, R. 2011, *ApJ*, 733, 32

* Gerola, H., & Seiden, P. E. 1978, *ApJ*, 223, 129

* Jokipii, J. R. 1966, *ApJ*, 146, 480

* Jokipii, J. R., & Lee, M. A. 2010, *ApJ*, 713, 475

* Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017, *ApJ*, 834, 25

* Kota, J., Merenyi, E., Jokipii, J. R., et al. 1982, *ApJ*, 254, 398

* Krause, M., Diehl, R., Böhringer, H., Freyberg, M., & Lubos, D. 2014, *A&A*, 566, A94

* Lerche, I., & Schlickeiser, R. 1985, *A&A*, 151, 408

* Melrose, D., & Crouch, A. 1997, *PASA*, 14, 251

* Parizot, E., Marcowith, A., van der Swaluw, E., Bykov, A. M., & Tatischeff, V. 2004, *A&A*, 424, 747

* Ptuskin, V. 1988, *SvAL*, 14, 255

* Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Heidelberg: Springer)

* Schlickeiser, R. 2009, *MPLA*, 24, 1461

* Schlickeiser, R. 2011, *ApJ*, 732, 96

* Schlickeiser, R. 2015, *PhPl*, 22, 091502

* Schneider, P. 1993, *A&A*, 278, 315

* Seo, J., Kang, H., & Ryu, D. 2018, *JKAS*, 51, 37

* Stecker, F. W. 1973, *ApJ*, 185, 499

* Ward, J. L., & Kruijssen, J. M. D. 2018, *MNRAS*, 475, 5659

* Weidner, C., Bonnell, I. A., & Zinnecker, H. 2010, *ApJ*, 724, 1503

* Wentzel, D. G. 1971, *ApJ*, 163, 503