regimes.

The time-dependent nature of our simulation means that protons escape the accelerator at different stages of the system evolution (i.e., as the SNR blast wave approaches the stellar wind termination shock) producing pions and neutrinos with varying hardness and maximum energy and these effects are included in our results. Of course, the complex evolution of the source and some unknown details of the mass distribution in the outer ISM region (e.g., the presence of dense shells or clouds) will also influence the results. However, we believe the general properties of our simulation are robust.

Assuming Bohm diffusion with $D_1(p) =D_2(p) = D(p) = cR_g(p)/3$, due to CR-driven, amplified magnetic instabilities in the CR accelerator, one obtains $$\label{TauAc1}
{t_\mathrm{acc}}\approx  \frac{c R_g(p)}{u_s\,u_w},$$ where $u_1=u_s$ is the SNR shock velocity, $u_2=u_w$ is the stellar wind speed, and $R_g(p)$ is the momentum dependent proton gyroradius. Then, using the scaling $B \approx \sqrt{4\pi \eta_b \rho}\,u_1$ for the amplified magnetic field, the acceleration time can be estimated as ${t_\mathrm{acc}}\approx 2\,\cdot10^{10}\,
\epsilon_{\rm PeV}\, (\eta_b n)^{-0.5}\,u_{s3}^{-2}\,u_{w3}^{-1}$ s. In the above expressions, $\rho=m_p n$, $\eta_b$ is the acceleration efficiency, i.e., the fraction of ram kinetic energy in the plasma flows converted into accelerated particles, the energy is in 